Index: trunk/psphot/src/psphotSourceSize.c
===================================================================
--- trunk/psphot/src/psphotSourceSize.c	(revision 27694)
+++ trunk/psphot/src/psphotSourceSize.c	(revision 27695)
@@ -50,9 +50,9 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-	if (i == chisqNum) continue; // skip chisq image
-	if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i, recipe, getPSFsize)) {
+        if (i == chisqNum) continue; // skip chisq image
+        if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i, recipe, getPSFsize)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on source size analysis for PSPHOT.INPUT entry %d", i);
-	    return false;
-	}
+            return false;
+        }
     }
     return true;
@@ -81,6 +81,6 @@
 
     if (!sources->n) {
-	psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
-	return true;
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
+        return true;
     }
 
@@ -140,5 +140,5 @@
     // XXX this should only be done on the first pass (ie, if we have newSources or allSources?)
     if (getPSFsize) {
-	psphotSourceSizePSF (&options, sources, psf);
+        psphotSourceSizePSF (&options, sources, psf);
     }
 
@@ -173,7 +173,9 @@
     pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
 
+    int num = 0;                        // Number of sources measured
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = sources->data[i];
         if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
+        num++;
 
         // replace object in image
@@ -200,4 +202,10 @@
         psVectorAppend (Ap, dMag);
         psVectorAppend (ApErr, source->errMag);
+    }
+    if (num == 0) {
+        // Not raising an error, because errors aren't being checked elsewhere in this function
+        psFree(Ap);
+        psFree(ApErr);
+        return false;
     }
 
@@ -258,5 +266,5 @@
             continue;
         }
-	// psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
+        // psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
     }
 
@@ -313,5 +321,5 @@
         }
 
-	// convert to Mmaj, Mmin:
+        // convert to Mmaj, Mmin:
         psF32 Mxx = source->moments->Mxx;
         psF32 Myy = source->moments->Myy;
@@ -338,26 +346,26 @@
         float dMag = source->psfMag - apMag;
 
-	// set nSigma to include both systematic and poisson error terms
-	// XXX the 'poisson error' contribution for size is probably wrong...
+        // set nSigma to include both systematic and poisson error terms
+        // XXX the 'poisson error' contribution for size is probably wrong...
         float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);
-	float nSigmaMXX = (Mxx - psfClump->X) / hypot(psfClump->dX, psfClump->X*psfClump->X*source->errMag);
-	float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag);
-
-	// partially-masked sources are more likely to be mis-measured PSFs
-	float sizeBias = 1.0;
-	if (source->pixWeight < 0.9) {
-	    sizeBias = 3.0;
-	} 
-
-	float minMxx = psfClump->X - sizeBias*options->nSigmaMoments*psfClump->dX;
-	float minMyy = psfClump->Y - sizeBias*options->nSigmaMoments*psfClump->dY;
-
-	// include MAG, MXX, and MYY?
+        float nSigmaMXX = (Mxx - psfClump->X) / hypot(psfClump->dX, psfClump->X*psfClump->X*source->errMag);
+        float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag);
+
+        // partially-masked sources are more likely to be mis-measured PSFs
+        float sizeBias = 1.0;
+        if (source->pixWeight < 0.9) {
+            sizeBias = 3.0;
+        }
+
+        float minMxx = psfClump->X - sizeBias*options->nSigmaMoments*psfClump->dX;
+        float minMyy = psfClump->Y - sizeBias*options->nSigmaMoments*psfClump->dY;
+
+        // include MAG, MXX, and MYY?
         source->extNsigma = nSigmaMAG;
 
-	// notes to clarify the source size classification rules:
-	// * a defect should be functionally equivalent to a cosmic ray
-	// * CR & defect should have a faintess limit (min S/N)
-	// * SAT stars should not be faint, but defects may?
+        // notes to clarify the source size classification rules:
+        // * a defect should be functionally equivalent to a cosmic ray
+        // * CR & defect should have a faintess limit (min S/N)
+        // * SAT stars should not be faint, but defects may?
 
         // Anything within this region is a probably PSF-like object. Saturated stars may land
@@ -365,10 +373,10 @@
         bool isPSF = (fabs(nSigmaMAG) < options->nSigmaApResid) && (fabs(nSigmaMXX) < sizeBias*options->nSigmaMoments) && (fabs(nSigmaMYY) < sizeBias*options->nSigmaMoments);
         if (isPSF) {
-	  psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
-		  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
-		  options->nSigmaApResid,sizeBias*options->nSigmaMoments);
-	  source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
-	  Npsf ++;
-	  continue;
+          psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
+                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
+                  options->nSigmaApResid,sizeBias*options->nSigmaMoments);
+          source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
+          Npsf ++;
+          continue;
         }
 
@@ -376,13 +384,13 @@
         // Defects may also be marked as SATSTAR -- XXX deactivate this flag?
         // XXX this rule is not great
-	// XXX only accept brightish detections as CRs
-	// (nSigmaMAG < -options->nSigmaApResid) || 
-	bool isCR = isCR = (source->errMag < 1.0 / SIZE_SN_LIM) && ((Mxx < minMxx) || (Myy < minMyy));
-	if (isCR) {
-	    psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g %f\t%g %g  %g %g  %g %g\t%g %g\t%g CR\t%g %g\n",
-		    source->peak->xf,source->peak->yf,source->pixWeight,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
-		    options->nSigmaApResid,sizeBias*options->nSigmaMoments);
+        // XXX only accept brightish detections as CRs
+        // (nSigmaMAG < -options->nSigmaApResid) ||
+        bool isCR = isCR = (source->errMag < 1.0 / SIZE_SN_LIM) && ((Mxx < minMxx) || (Myy < minMyy));
+        if (isCR) {
+            psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g %f\t%g %g  %g %g  %g %g\t%g %g\t%g CR\t%g %g\n",
+                    source->peak->xf,source->peak->yf,source->pixWeight,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
+                    options->nSigmaApResid,sizeBias*options->nSigmaMoments);
             source->mode |= PM_SOURCE_MODE_DEFECT;
-	    source->tmpFlags |= PM_SOURCE_TMPF_SIZE_CR_CANDIDATE;
+            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_CR_CANDIDATE;
             Ncr ++;
             continue;
@@ -392,8 +400,8 @@
         // just large saturated regions.
         if (source->mode & PM_SOURCE_MODE_SATSTAR) {
-	    psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g SAT\t%g %g\n",
-		    source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
-		    options->nSigmaApResid,sizeBias*options->nSigmaMoments);
-	    source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
+            psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g SAT\t%g %g\n",
+                    source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
+                    options->nSigmaApResid,sizeBias*options->nSigmaMoments);
+            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
             Nsat ++;
             continue;
@@ -403,20 +411,20 @@
         bool isEXT = (nSigmaMAG > options->nSigmaApResid) || (Mxx > minMxx) || (Myy > minMyy);
         if (isEXT) {
-	  psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Ext\t%g %g\n",
-		  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
-		  options->nSigmaApResid,sizeBias*options->nSigmaMoments);
-	  
+          psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Ext\t%g %g\n",
+                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
+                  options->nSigmaApResid,sizeBias*options->nSigmaMoments);
+
             source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
-	    source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
+            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
             Next ++;
             continue;
         }
-	psTrace("psphotSourceClassRegion.MISS",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Unk\t%g %g\n",
-		source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
-		options->nSigmaApResid,sizeBias*options->nSigmaMoments);
-	
-	// sources that reach here are probably too faint for a reasonable source size measurement
+        psTrace("psphotSourceClassRegion.MISS",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Unk\t%g %g\n",
+                source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
+                options->nSigmaApResid,sizeBias*options->nSigmaMoments);
+
+        // sources that reach here are probably too faint for a reasonable source size measurement
         // psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigmaMAG);
-	source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
+        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
         Nmiss ++;
     }
@@ -450,6 +458,6 @@
 
         // skip unless this source is thought to be a cosmic ray.  flag the detection and mask the pixels
-	// XXX this may be degenerate with the above test
-	if (!(source->mode & PM_SOURCE_MODE_DEFECT)) continue;
+        // XXX this may be degenerate with the above test
+        if (!(source->mode & PM_SOURCE_MODE_DEFECT)) continue;
 
         // Integer position of peak
@@ -463,8 +471,8 @@
             continue;
         }
-	
-	// XXX for testing, only CRMASK a single source:
-	if (options->xtest && (fabs(source->peak->xf - options->xtest) > 5)) continue;
-	if (options->ytest && (fabs(source->peak->yf - options->ytest) > 5)) continue;
+
+        // XXX for testing, only CRMASK a single source:
+        if (options->xtest && (fabs(source->peak->xf - options->xtest) > 5)) continue;
+        if (options->ytest && (fabs(source->peak->yf - options->ytest) > 5)) continue;
 
         // replace object in image
@@ -473,17 +481,17 @@
         }
 
-	// XXX this is running slowly and is too agressive, but it more-or-less works
-	psTrace("psphot", 6, "mask cosmic ray at %f, %f\n", source->peak->xf, source->peak->yf);
-	if (options->apply) {
-	    psphotMaskCosmicRay(readout, source, options->crMask);
-	} else {
-	    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
-	}
-	nMasked ++;
+        // XXX this is running slowly and is too agressive, but it more-or-less works
+        psTrace("psphot", 6, "mask cosmic ray at %f, %f\n", source->peak->xf, source->peak->yf);
+        if (options->apply) {
+            psphotMaskCosmicRay(readout, source, options->crMask);
+        } else {
+            source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+        }
+        nMasked ++;
 
         // re-subtract the object, leave local sky
         pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
     }
-    
+
     // now that we have masked pixels associated with CRs, we can grow the mask
     if (options->grow > 0) {
@@ -503,5 +511,5 @@
     // XXX test : save the mask image
     if (0) {
-	psphotSaveImage (NULL, readout->mask,   "mask.fits");
+        psphotSaveImage (NULL, readout->mask,   "mask.fits");
     }
 
@@ -540,5 +548,5 @@
     psImage *image= readout->image;
     psImage *variance = readout->variance;
-      
+
     int binning = 2;
     float sigma_thresh = 3.0;
@@ -553,104 +561,104 @@
     psImage *edges  = psImageAlloc(dx,dy,image->type.type);
     psImage *mymask = psImageAlloc(dx,dy,PS_TYPE_IMAGE_MASK);
-      
+
     // Load my copy of things.
     for (int y = 0; y < dy; y++) {
-	for (int x = 0; x < dx; x++) {
-	    mypix->data.F32[y][x] = image->data.F32[y+ys][x+xs];
-	    myvar->data.F32[y][x] = variance->data.F32[y+ys][x+xs];
-	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0x00;
-	}
+        for (int x = 0; x < dx; x++) {
+            mypix->data.F32[y][x] = image->data.F32[y+ys][x+xs];
+            myvar->data.F32[y][x] = variance->data.F32[y+ys][x+xs];
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0x00;
+        }
     }
     // Mask so I can see on the output image where the footprint is.
     for (int i = 0; i < footprint->spans->n; i++) {
-	pmSpan *sp = footprint->spans->data[i];
-	for (int j = sp->x0; j <= sp->x1; j++) {
-	    int y = sp->y - ys;
-	    int x = j - xs;
-	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x01;
-	}
+        pmSpan *sp = footprint->spans->data[i];
+        for (int j = sp->x0; j <= sp->x1; j++) {
+            int y = sp->y - ys;
+            int x = j - xs;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x01;
+        }
     }
 
     int nCRpix = 1; // force at least one pass...
     for (int iteration = 0; (iteration < max_iter) && (nCRpix > 0); iteration++) {
-	nCRpix = 0;
-	psImageInit (binned, 0.0);
-	psImageInit (conved, 0.0);
-	psImageInit (edges, 0.0);
-
-	// Make subsampled image. Maybe this should be called "unbinned" or something
-	for (int y = 0; y < binning * dy; y++) {
-	    int yraw = y / binning;
-	    for (int x = 0; x < binning * dx; x++) {
-		int xraw = x / binning;
-		binned->data.F32[y][x] = mypix->data.F32[yraw][xraw];
-	    }
-	}
-
-	// Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero
-	for (int y = 1; y < binning * dy - 1; y++) {
-	    for (int x = 1; x < binning * dx - 1; x++) {
-		float value = binned->data.F32[y][x] - 0.25 * 
-		    (binned->data.F32[y+0][x-1] + binned->data.F32[y+0][x+1] + 
-		     binned->data.F32[y-1][x+0] + binned->data.F32[y+1][x+0]);
-		value = PS_MAX(0.0, value);
-
-		conved->data.F32[y][x] = value;
-	    }
-	}
-
-	// Create an edge map by rebinning
-	for (int y = 0; y < binning * dy; y++) {
-	    int yraw = y / binning;
-	    for (int x = 0; x < binning * dx; x++) {
-		int xraw = x / binning;
-		edges->data.F32[yraw][xraw] += conved->data.F32[y][x];
-	    }
-	}
-
-	// coordinate of peak in subimage pixels:
-	int xPeak = peak->x - xs;
-	int yPeak = peak->y - ys;
-
-	// Modify my mask if we're above the significance threshold, but only for connected pixels
-	nCRpix = psphotMaskCosmicRayConnected (xPeak, yPeak, mymask, myvar, edges, binning, sigma_thresh);
-	
+        nCRpix = 0;
+        psImageInit (binned, 0.0);
+        psImageInit (conved, 0.0);
+        psImageInit (edges, 0.0);
+
+        // Make subsampled image. Maybe this should be called "unbinned" or something
+        for (int y = 0; y < binning * dy; y++) {
+            int yraw = y / binning;
+            for (int x = 0; x < binning * dx; x++) {
+                int xraw = x / binning;
+                binned->data.F32[y][x] = mypix->data.F32[yraw][xraw];
+            }
+        }
+
+        // Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero
+        for (int y = 1; y < binning * dy - 1; y++) {
+            for (int x = 1; x < binning * dx - 1; x++) {
+                float value = binned->data.F32[y][x] - 0.25 *
+                    (binned->data.F32[y+0][x-1] + binned->data.F32[y+0][x+1] +
+                     binned->data.F32[y-1][x+0] + binned->data.F32[y+1][x+0]);
+                value = PS_MAX(0.0, value);
+
+                conved->data.F32[y][x] = value;
+            }
+        }
+
+        // Create an edge map by rebinning
+        for (int y = 0; y < binning * dy; y++) {
+            int yraw = y / binning;
+            for (int x = 0; x < binning * dx; x++) {
+                int xraw = x / binning;
+                edges->data.F32[yraw][xraw] += conved->data.F32[y][x];
+            }
+        }
+
+        // coordinate of peak in subimage pixels:
+        int xPeak = peak->x - xs;
+        int yPeak = peak->y - ys;
+
+        // Modify my mask if we're above the significance threshold, but only for connected pixels
+        nCRpix = psphotMaskCosmicRayConnected (xPeak, yPeak, mymask, myvar, edges, binning, sigma_thresh);
+
 # if DUMPPICS
-	psphotSaveImage (NULL, mypix,   "crmask.pix.fits");
+        psphotSaveImage (NULL, mypix,   "crmask.pix.fits");
 # endif
-	
+
 // XXX do not repair the pixels in isophot version
 # if 0
-	// "Repair" Masked pixels for the next round.
-	for (int y = 1; y < dy - 1; y++) {
-	    for (int x = 1; x < dx - 1; x++) {
-		if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) {
-		    myfix->data.F32[y][x] = mypix->data.F32[y][x];
-		    continue;
-		}
-		myfix->data.F32[y][x] = 0.25 * 
-		    (mypix->data.F32[y+0][x-1] + mypix->data.F32[y+0][x+1] +
-		     mypix->data.F32[y-1][x+0] + mypix->data.F32[y+1][x+0]);
-	    }
-	}
-	
-	// "Repair" Masked pixels for the next round.
-	for (int y = 1; y < dy - 1; y++) {
-	    for (int x = 1; x < dx - 1; x++) {
-		mypix->data.F32[y][x] = myfix->data.F32[y][x];
-	    }
-	}
+        // "Repair" Masked pixels for the next round.
+        for (int y = 1; y < dy - 1; y++) {
+            for (int x = 1; x < dx - 1; x++) {
+                if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) {
+                    myfix->data.F32[y][x] = mypix->data.F32[y][x];
+                    continue;
+                }
+                myfix->data.F32[y][x] = 0.25 *
+                    (mypix->data.F32[y+0][x-1] + mypix->data.F32[y+0][x+1] +
+                     mypix->data.F32[y-1][x+0] + mypix->data.F32[y+1][x+0]);
+            }
+        }
+
+        // "Repair" Masked pixels for the next round.
+        for (int y = 1; y < dy - 1; y++) {
+            for (int x = 1; x < dx - 1; x++) {
+                mypix->data.F32[y][x] = myfix->data.F32[y][x];
+            }
+        }
 # endif
 
 # if DUMPPICS
-	fprintf (stderr, "CRMASK %d %d %d %d %d\n", xs, ys, dx, dy, iteration);
-	psphotSaveImage (NULL, mypix,   "crmask.fix.fits");
-	psphotSaveImage (NULL, myvar,   "crmask.var.fits");
-	psphotSaveImage (NULL, binned,  "crmask.binn.fits");
-	psphotSaveImage (NULL, conved,  "crmask.conv.fits");
-	psphotSaveImage (NULL, edges,   "crmask.edge.fits");
-	psphotSaveImage (NULL, mymask,  "crmask.mask.fits");
+        fprintf (stderr, "CRMASK %d %d %d %d %d\n", xs, ys, dx, dy, iteration);
+        psphotSaveImage (NULL, mypix,   "crmask.fix.fits");
+        psphotSaveImage (NULL, myvar,   "crmask.var.fits");
+        psphotSaveImage (NULL, binned,  "crmask.binn.fits");
+        psphotSaveImage (NULL, conved,  "crmask.conv.fits");
+        psphotSaveImage (NULL, edges,   "crmask.edge.fits");
+        psphotSaveImage (NULL, mymask,  "crmask.mask.fits");
 # endif
-	psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration, nCRpix);
+        psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration, nCRpix);
     }
 
@@ -659,20 +667,20 @@
     // XXX can't we use nCRpix == 1 to test for these?
     for (int x = 0; x < dx; x++) {
-	for (int y = 0; y < dy; y++) {
-	    if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) continue;
-	    if ((x-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) {
-		continue;
-	    }
-	    if ((y-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) {
-		continue;
-	    }
-	    if ((x+1 < dx) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) {
-		continue;
-	    }
-	    if ((y+1 < dy) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) {
-		continue;
-	    }
-	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40;
-	}
+        for (int y = 0; y < dy; y++) {
+            if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) continue;
+            if ((x-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) {
+                continue;
+            }
+            if ((y-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) {
+                continue;
+            }
+            if ((x+1 < dx) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) {
+                continue;
+            }
+            if ((y+1 < dy) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) {
+                continue;
+            }
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40;
+        }
     }
 # endif
@@ -681,16 +689,16 @@
     nCRpix = 0;
     for (int x = 0; x < dx; x++) {
-	for (int y = 0; y < dy; y++) {
-	    if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
-		mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ys+mask->row0][x+xs+mask->col0] |= maskVal;
-		nCRpix ++;
-	    }
-	}
-    }
-      
+        for (int y = 0; y < dy; y++) {
+            if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
+                mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ys+mask->row0][x+xs+mask->col0] |= maskVal;
+                nCRpix ++;
+            }
+        }
+    }
+
     // XXX if we decide this REALLY is a cosmic ray, set the CR_LIMIT bit
     if (nCRpix > 1) {
-	source->mode |= PM_SOURCE_MODE_CR_LIMIT;
-	source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
+        source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     }
     // fprintf (stderr, "CRMASK %d %d %d %d %d\n", peak->x, peak->y, dx, dy, nCRpix);
@@ -703,5 +711,5 @@
     psFree(edges);
     psFree(mymask);
-      
+
     return true;
 }
@@ -711,11 +719,11 @@
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = sources->data[i];
-	pmPeak *peak = source->peak;
-	pmFootprint *footprint = peak->footprint;
-	if (!footprint) continue;
-	for (int j = 0; j < footprint->spans->n; j++) {
-	    pmSpan *sp = footprint->spans->data[j];
-	    psAssert (sp, "missing span");
-	}
+        pmPeak *peak = source->peak;
+        pmFootprint *footprint = peak->footprint;
+        if (!footprint) continue;
+        for (int j = 0; j < footprint->spans->n; j++) {
+            pmSpan *sp = footprint->spans->data[j];
+            psAssert (sp, "missing span");
+        }
     }
     return true;
@@ -741,5 +749,5 @@
     if (!footprint) {
       psTrace("psphot.czw",2,"Using isophot CR mask code.");
-      
+
         // if we have not footprint, use the old code to mask by isophot
         psphotMaskCosmicRayIsophot (source, maskVal, crMask);
@@ -749,5 +757,5 @@
     if (!footprint->spans) {
       psTrace("psphot.czw",2,"Using isophot CR mask code.");
-      
+
         // if we have no footprint, use the old code to mask by isophot
         psphotMaskCosmicRayIsophot (source, maskVal, crMask);
@@ -756,15 +764,15 @@
     psphotMaskCosmicRayIsophot (source, maskVal, crMask);
     // mask all of the pixels covered by the spans of the footprint
-    for (int j = 1; j < footprint->spans->n; j++) { 
-	pmSpan *span1 = footprint->spans->data[j]; 
-
-	int iy = span1->y; 
-	int xs = span1->x0; 
-	int xe = span1->x1; 
-
-	for (int ix = xs; ix < xe; ix++) { 
-	    mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; 
-	} 
-    } 
+    for (int j = 1; j < footprint->spans->n; j++) {
+        pmSpan *span1 = footprint->spans->data[j];
+
+        int iy = span1->y;
+        int xs = span1->x0;
+        int xe = span1->x1;
+
+        for (int ix = xs; ix < xe; ix++) {
+            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
+        }
+    }
     return true;
 }
@@ -772,34 +780,34 @@
 # define VERBOSE 0
 int psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh) {
-    
+
     int xLo, xRo;
     int nCRpix = 0;
 
     float noise_factor = 5.0 / 4.0;  // Intrinsic to the Laplacian making noise spikes spikier.
-    
+
     // mark the pixels in this row to the left, then the right. stay within footprint
     int xL = xPeak; // find the range of valid pixels in this row
     int xR = xPeak;
     for (int ix = xPeak; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] & 0x01); ix--) {
-	float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]);
-	float value = edges->data.F32[yPeak][ix] / noise;
-	if (value < sigma_thresh ) break;
-	mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40;
-	xL = ix;
-	nCRpix ++;
-	if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR);
+        float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]);
+        float value = edges->data.F32[yPeak][ix] / noise;
+        if (value < sigma_thresh ) break;
+        mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40;
+        xL = ix;
+        nCRpix ++;
+        if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR);
     }
     for (int ix = xPeak; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] & 0x01); ix++) {
-	float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]);
-	float value = edges->data.F32[yPeak][ix] / noise;
-	if (value < sigma_thresh ) break;
-	mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40;
-	xR = ix;
-	nCRpix ++;
-	if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR);
+        float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]);
+        float value = edges->data.F32[yPeak][ix] / noise;
+        if (value < sigma_thresh ) break;
+        mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40;
+        xR = ix;
+        nCRpix ++;
+        if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR);
     }
     // xL and xR mark the first and last valid pixel in the row
 
-    // for each of the neighboring rows, mark the high pixels if they touch the range xL to xR 
+    // for each of the neighboring rows, mark the high pixels if they touch the range xL to xR
     xLo = PS_MAX(xL - 1, 0);
     xRo = PS_MIN(xR + 1, mymask->numCols);
@@ -808,50 +816,50 @@
     for (int iy = yPeak - 1; iy >= 0; iy--) {
 
-	int xLn = -1;
-	int xRn = -1;
-	int newPix = 0;
-
-	// mark the pixels in the good range
-	for (int ix = xLo; ix < xRo; ix++) {
-	    if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint
- 	    float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
-	    float value = edges->data.F32[iy][ix] / noise;
-	    if (value < sigma_thresh ) continue;
-	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
-	    if (xLn == -1) xLn = ix; // first valid pixel in this row
-	    xRn = ix;		     // last valid pixel in this row
-	    nCRpix ++;
-	    newPix ++;
-	    if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
-	}
-
-	// mark the pixels to the left of the good range
-	for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) {
-	    float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
-	    float value = edges->data.F32[iy][ix] / noise;
-	    if (value < sigma_thresh ) break;
-	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
-	    if (xRn == -1) xRn = ix; // last valid pixel in this row
-	    xLn = ix;
-	    nCRpix ++;
-	    newPix ++;
-	    if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
-	}
-
-	// mark the pixels to the right of the good range
-	for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) {
-	    float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
-	    float value = edges->data.F32[iy][ix] / noise;
-	    if (value < sigma_thresh ) break;
-	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
-	    if (xLn == -1) xLn = ix; // first valid pixel in this row
-	    xRn = ix;
-	    nCRpix ++;
-	    newPix ++;
-	    if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
-	}
-	if (newPix == 0) break;
-	xLo = PS_MAX(xLn - 1, 0);
-	xRo = PS_MIN(xRn + 1, mymask->numCols);
+        int xLn = -1;
+        int xRn = -1;
+        int newPix = 0;
+
+        // mark the pixels in the good range
+        for (int ix = xLo; ix < xRo; ix++) {
+            if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) continue;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xLn == -1) xLn = ix; // first valid pixel in this row
+            xRn = ix;                // last valid pixel in this row
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+
+        // mark the pixels to the left of the good range
+        for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) {
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) break;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xRn == -1) xRn = ix; // last valid pixel in this row
+            xLn = ix;
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+
+        // mark the pixels to the right of the good range
+        for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) {
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) break;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xLn == -1) xLn = ix; // first valid pixel in this row
+            xRn = ix;
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+        if (newPix == 0) break;
+        xLo = PS_MAX(xLn - 1, 0);
+        xRo = PS_MIN(xRn + 1, mymask->numCols);
     }
 
@@ -862,50 +870,50 @@
     for (int iy = yPeak + 1; iy < mymask->numRows; iy++) {
 
-	int xLn = -1;
-	int xRn = -1;
-	int newPix = 0;
-
-	// mark the pixels in the good range
-	for (int ix = xLo; ix < xRo; ix++) {
-	    if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint
- 	    float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
-	    float value = edges->data.F32[iy][ix] / noise;
-	    if (value < sigma_thresh ) continue;
-	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
-	    if (xLn == -1) xLn = ix; // first valid pixel in this row
-	    xRn = ix;		     // last valid pixel in this row
-	    nCRpix ++;
-	    newPix ++;
-	    if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
-	}
-
-	// mark the pixels to the left of the good range
-	for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) {
-	    float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
-	    float value = edges->data.F32[iy][ix] / noise;
-	    if (value < sigma_thresh ) break;
-	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
-	    if (xRn == -1) xRn = ix; // last valid pixel in this row
-	    xLn = ix;
-	    nCRpix ++;
-	    newPix ++;
-	    if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
-	}
-
-	// mark the pixels to the right of the good range
-	for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) {
-	    float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
-	    float value = edges->data.F32[iy][ix] / noise;
-	    if (value < sigma_thresh ) break;
-	    mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
-	    if (xLn == -1) xLn = ix; // first valid pixel in this row
-	    xRn = ix;
-	    nCRpix ++;
-	    newPix ++;
-	    if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
-	}
-	if (newPix == 0) break;
-	xLo = PS_MAX(xLn - 1, 0);
-	xRo = PS_MIN(xRn + 1, mymask->numCols);
+        int xLn = -1;
+        int xRn = -1;
+        int newPix = 0;
+
+        // mark the pixels in the good range
+        for (int ix = xLo; ix < xRo; ix++) {
+            if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) continue;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xLn == -1) xLn = ix; // first valid pixel in this row
+            xRn = ix;                // last valid pixel in this row
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+
+        // mark the pixels to the left of the good range
+        for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) {
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) break;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xRn == -1) xRn = ix; // last valid pixel in this row
+            xLn = ix;
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+
+        // mark the pixels to the right of the good range
+        for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) {
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) break;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xLn == -1) xLn = ix; // first valid pixel in this row
+            xRn = ix;
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+        if (newPix == 0) break;
+        xLo = PS_MAX(xLn - 1, 0);
+        xRo = PS_MIN(xRn + 1, mymask->numCols);
     }