Index: /trunk/psphot/.cvsignore
===================================================================
--- /trunk/psphot/.cvsignore	(revision 4129)
+++ /trunk/psphot/.cvsignore	(revision 4129)
@@ -0,0 +1,1 @@
+bin
Index: /trunk/psphot/Makefile
===================================================================
--- /trunk/psphot/Makefile	(revision 4128)
+++ /trunk/psphot/Makefile	(revision 4129)
@@ -32,6 +32,10 @@
 $(SRC)/mark_psf_sources.$(ARCH).o \
 $(SRC)/subtract_psf_sources.$(ARCH).o \
+$(SRC)/mark_psf_source.$(ARCH).o \
+$(SRC)/subtract_psf_source.$(ARCH).o \
+$(SRC)/by_SN.$(ARCH).o \
 $(SRC)/fit_galaxies.$(ARCH).o \
-$(SRC)/subtract_galaxies.$(ARCH).o
+$(SRC)/subtract_galaxies.$(ARCH).o \
+$(SRC)/find_defects.$(ARCH).o
 
 psphot: $(BIN)/psphot.$(ARCH)
Index: /trunk/psphot/src/apply_psf_model.c
===================================================================
--- /trunk/psphot/src/apply_psf_model.c	(revision 4128)
+++ /trunk/psphot/src/apply_psf_model.c	(revision 4129)
@@ -25,5 +25,5 @@
 
     float shapeNsigma = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA");
-    float snFaint     = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
+    // float snFaint     = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
 
     // set the object surface-brightness limit for fitted pixels
@@ -66,5 +66,5 @@
 	Nfit ++;
 
-	mark_psf_source (source);
+	mark_psf_source (source, shapeNsigma);
 	subtract_psf_source (source);
     }
Index: /trunk/psphot/src/by_SN.c
===================================================================
--- /trunk/psphot/src/by_SN.c	(revision 4128)
+++ /trunk/psphot/src/by_SN.c	(revision 4129)
@@ -1,8 +1,9 @@
 # include "psphot.h"
 
+// sort by SN (descending)
 int by_SN (const void **a, const void **b)
 {
-    psSource *A = *a;
-    psSource *B = *b;
+    psSource *A = *(psSource **)a;
+    psSource *B = *(psSource **)b;
 
     psF32 fA = (A->moments == NULL) ? 0 : A->moments->SN;
@@ -10,6 +11,6 @@
 
     psF32 diff = fA - fB;
-    if (diff < FLT_EPSILON) return (-1);
-    if (diff > FLT_EPSILON) return (+1);
+    if (diff > FLT_EPSILON) return (-1);
+    if (diff < FLT_EPSILON) return (+1);
     return (0);
 }
Index: /trunk/psphot/src/choose_psf_model.c
===================================================================
--- /trunk/psphot/src/choose_psf_model.c	(revision 4128)
+++ /trunk/psphot/src/choose_psf_model.c	(revision 4129)
@@ -20,4 +20,5 @@
 	psArrayAdd (stars, 200, source);
     }
+    psTrace (".psphot.pspsf", 3, "selected candidate %d PSF objects\n", stars->n);
 
     // define model fit pixels
Index: /trunk/psphot/src/find_defects.c
===================================================================
--- /trunk/psphot/src/find_defects.c	(revision 4129)
+++ /trunk/psphot/src/find_defects.c	(revision 4129)
@@ -0,0 +1,121 @@
+# include "psphot.h"
+# define NSIGMA 3.0
+
+void dump_image (psImage *image, char *filename);
+
+bool find_defects (psImage *zapmask, psArray *sources, psMetadata *config, pmPSF *psf) 
+{ 
+
+  bool status;
+  float gradx, grady, gradlim;
+  float N0, N1, peak, sky;
+  int x, y;
+
+  // measure the x and y max gradients
+  float RADIUS = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
+  int dN = 2*RADIUS + 1;
+
+  psImage *flux = psImageAlloc (dN, dN, PS_TYPE_F32);
+
+  // create an object located at the image center
+  psModel *model = psModelAlloc (psf->type);
+  model->params->data.F32[0] = 0;
+  model->params->data.F32[1] = 1;
+  // XXX : temporary hack : center of 2k x 4k image
+  model->params->data.F32[2] = 1024;
+  model->params->data.F32[3] = 2048;
+
+  // fill in the PSF parameters
+  psModelFromPSFFunc modelFromPSFFunc = psModelFromPSFFunc_GetFunction (psf->type);
+  modelFromPSFFunc (model, model, psf);
+
+  // generate the image flux
+  pmSourceAddModel (flux, NULL, model, true);
+  dump_image (flux, "psf.fits");
+
+  fprintf (stderr, "sx: %f, sy: %f\n", model->params->data.F32[4], model->params->data.F32[5]);
+
+  // measure max dfdx
+  gradx = 0;
+  for (int iy = 0; iy < flux->numRows; iy++) {
+    for (int ix = 0; ix < flux->numCols - 1; ix++) {
+      gradx = PS_MAX (gradx, fabs(flux->data.F32[iy][ix] - flux->data.F32[iy][ix + 1]));
+    }
+  }      
+
+  // measure max dfdy
+  grady = 0;
+  for (int iy = 0; iy < flux->numRows - 1; iy++) {
+    for (int ix = 0; ix < flux->numCols; ix++) {
+      grady = PS_MAX (grady, fabs(flux->data.F32[iy][ix] - flux->data.F32[iy + 1][ix]));
+    }
+  }      
+  psLogMsg ("psphot.defects", 4, "gradx: %f, grady: %f\n", gradx, grady);
+  psFree (flux);
+  psFree (model);
+      
+  // search for pixels with (dfdy > grady) or (dfdx > gradx) allowing for Nsigma outliers
+  for (int i = 0; i < sources->n; i++) {
+    psSource *source = sources->data[i];
+    psImage *image = source->pixels;
+    psImage *mask  = source->mask;
+    psImage *zap = psImageAlloc (image->numCols, image->numRows, PS_TYPE_U8);
+    psImageInit (zap, 0);
+
+    // get sky-subtracted peak (skip sources with insufficient gradient)
+    sky  = source->moments->Sky;
+    peak = source->moments->Peak - sky;
+    if (peak < peak*PS_MIN(gradx, grady) + NSIGMA*sqrt(2*sky + peak)) continue;
+    psTrace (".psphot.defects", 4, "peak: %f, SN: %f\n", peak, source->moments->SN);
+
+    // search for (dfdx)
+    for (int iy = 0; iy < image->numRows; iy++) {
+      for (int ix = 0; ix < image->numCols - 1; ix++) {
+	N0 = image->data.F32[iy][ix];
+	N1 = image->data.F32[iy][ix + 1];
+	gradlim = fabs(N1 - N0 - NSIGMA*(sqrt(N1) + sqrt(N0))) / peak;
+	// XXX this line is wrong: we should use the appropriate noise from the noise map
+	
+	// if gradient is too large, zap the high pixel
+	if (gradlim < gradx) continue;
+	if (N1 > N0) {
+	  zap->data.U8[iy][ix + 1] = 1;
+	} else {
+	  zap->data.U8[iy][ix] = 1;
+	}
+      }
+    }
+
+    // search for (dfdy)
+    for (int iy = 0; iy < image->numRows - 1; iy++) {
+      for (int ix = 0; ix < image->numCols; ix++) {
+	N0 = image->data.F32[iy][ix];
+	N1 = image->data.F32[iy + 1][ix];
+	gradlim = fabs(N1 - N0 - NSIGMA*(sqrt(N1) + sqrt(N0))) / peak;
+	// XXX this line is wrong: we should use the appropriate noise from the noise map
+	
+	// if gradient is too large, zap the high pixel
+	// note the inefficiency: we are re-testing known zapped entries
+	if (gradlim > grady) continue;
+	if (N1 > N0) {
+	  zap->data.U8[iy + 1][ix] = 1;
+	} else {
+	  zap->data.U8[iy][ix] = 1;
+	}	    
+      }
+    }
+    psFree (zap);
+
+    // mask & zero the zapped pixels
+    for (int iy = 0; iy < image->numRows - 1; iy++) {
+      for (int ix = 0; ix < image->numCols; ix++) {
+	if (!zap->data.U8[iy][ix]) continue;
+	x = ix + image->col0;
+	y = iy + image->row0;
+	zapmask->data.F32[y][x] = 1;
+      }
+    }
+
+  }
+  return (true);
+}
Index: /trunk/psphot/src/mark_psf_source.c
===================================================================
--- /trunk/psphot/src/mark_psf_source.c	(revision 4128)
+++ /trunk/psphot/src/mark_psf_source.c	(revision 4129)
@@ -14,7 +14,6 @@
 // any object which meets the criterion is marked as 
 // PS_SOURCE_BRIGHTSTAR 
-bool mark_psf_source (psArray *sources)
+bool mark_psf_source (psSource *source, float shapeNsigma)
 { 
-    bool  status      = false;
     float dSX, dSY, SX, SY, SN;
     float nSx, nSy;
Index: /trunk/psphot/src/psPolynomials.c
===================================================================
--- /trunk/psphot/src/psPolynomials.c	(revision 4128)
+++ /trunk/psphot/src/psPolynomials.c	(revision 4129)
@@ -233,5 +233,5 @@
     // Build the B and A data structs.
     for (int k = 0; k < x->n; k++) {
-        if ((mask != NULL) && mask->U8[k]) continue;
+        if ((mask != NULL) && mask->data.U8[k]) continue;
 	xSums = BuildSums1D(xSums, x->data.F64[k], nTerm);
       
@@ -325,5 +325,5 @@
     // Build the B and A data structs.
     for (int k  = 0; k < x->n; k++) {
-      if ((mask != NULL) && mask->data.U8[i]) continue;
+      if ((mask != NULL) && mask->data.U8[k]) continue;
 	Sums = BuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm);
       
Index: /trunk/psphot/src/psUtils.c
===================================================================
--- /trunk/psphot/src/psUtils.c	(revision 4128)
+++ /trunk/psphot/src/psUtils.c	(revision 4129)
@@ -261,2 +261,50 @@
 }
     
+bool psImageInit (psImage *image,...) {
+
+  va_list argp;
+  psU8  vU8;
+  psF32 vF32;
+  psF64 vF64;
+
+  if (image == NULL) return (false);
+
+  va_start (argp, image);
+
+  switch (image->type.type) {
+    case PS_TYPE_U8:
+      vU8 = va_arg (argp, psU32);
+
+      for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	  image->data.U8[iy][ix] = vU8;
+	}
+      }
+      break;
+
+    case PS_TYPE_F32:
+      vF32 = va_arg (argp, psF64);
+      
+      for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	  image->data.F32[iy][ix] = vF32;
+	}
+      }
+      return (true);
+
+    case PS_TYPE_F64:
+      vF64 = va_arg (argp, psF64);
+
+      for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	  image->data.F64[iy][ix] = vF64;
+	}
+      }
+      return (true);
+
+    default:
+      psAbort ("psphot.psUtils", "datatype %d not defined in psImageInit\n", image->type);
+      return (false);
+  }
+  return (false);
+}	    
Index: /trunk/psphot/src/psphot.c
===================================================================
--- /trunk/psphot/src/psphot.c	(revision 4128)
+++ /trunk/psphot/src/psphot.c	(revision 4129)
@@ -1,3 +1,4 @@
 # include "psphot.h"
+void dump_image (psImage *image, char *filename);
 
 int main (int argc, char **argv) {
@@ -33,4 +34,10 @@
     psf = choose_psf_model (image, config, sources);
 
+    psImage *zap = psImageAlloc (image->numCols, image->numRows, PS_TYPE_F32);
+    psImageInit (zap, 0);
+    find_defects (zap, sources, config, psf);
+    dump_image (zap, argv[2]);
+    exit (0);
+
     // test PSF on all except SATURATE and DEFECT (artifacts)
     apply_psf_model (image, config, sources, psf, sky);
@@ -44,10 +51,5 @@
     // subtract_galaxies (sources, config);
  
-    // write out subtracted image
-    unlink (argv[2]);
-    psFits *fits = psFitsAlloc (argv[2]);
-    psFitsWriteImage (fits, NULL, image, 0, NULL);
-    psFree (fits);
-
+    dump_image (image, argv[2]);
     DumpModelPSF (sources, "sources.dat");
     exit (0);
@@ -59,2 +61,12 @@
     exit (2);
 }
+
+void dump_image (psImage *image, char *filename) {
+
+    unlink (filename);
+    psFits *fits = psFitsAlloc (filename);
+    psFitsWriteImage (fits, NULL, image, 0, NULL);
+    psFree (fits);
+    return;
+}
+
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 4128)
+++ /trunk/psphot/src/psphot.h	(revision 4129)
@@ -21,4 +21,5 @@
 psF64 psTimerMark (char *name);
 
+bool psImageInit (psImage *image,...);
 void psImageSmooth (psImage *image, float sigma, float Nsigma);
 psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in);
@@ -40,9 +41,10 @@
 psF32 Polynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y);
 psImage *psBuildSums2D(psImage* sums,psF64 x,psF64 y,psS32 nXterm, psS32 nYterm);
-psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, const psVector* x, const psVector* y, const psVector* yErr);
-psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
+psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* yErr);
+psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
 psPolynomial2D* Polynomial2DAlloc(psS32 nXorder, psS32 nYorder, psPolynomialType type);
 void psPolynomial2DDump (psPolynomial2D *poly);
-psPolynomial2D* RobustFit2D(psPolynomial2D* poly, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
+psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
+psPolynomial2D* RobustFit2D(psPolynomial2D* poly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
 psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly, const psVector *x,const psVector *y);
 psVector *psBuildSums1D(psVector* sums, psF64 x,psS32 nTerm);
@@ -58,5 +60,5 @@
 pmPSF_Test *pmPSF_TestModel (psArray *stars, char *modelName);
 
-bool pmPSFFromModels (pmPSF *psf, psArray *models);
+bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask);
 psModel *psModelFromPSF (psModel *model, pmPSF *psf);
 void pmSourceMaskRegion (psSource *source, psRegion *region);
@@ -75,2 +77,7 @@
 bool fit_galaxies (psImage *image, psMetadata *config, psArray *sources);
 bool subtract_galaxies (psArray *sources, psMetadata *config);
+
+int by_SN (const void **a, const void **b);
+bool subtract_psf_source (psSource *source);
+bool mark_psf_source (psSource *source, float shapeNsigma);
+bool find_defects (psImage *zapmask, psArray *sources, psMetadata *config, pmPSF *psf);
Index: /trunk/psphot/src/pspsf.c
===================================================================
--- /trunk/psphot/src/pspsf.c	(revision 4128)
+++ /trunk/psphot/src/pspsf.c	(revision 4129)
@@ -68,8 +68,17 @@
     }
     psLogMsg ("psphot.psftest", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
+    psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
     DumpModelFits (test->modelFLT, "modelsFLT.dat");
 
     // stage 2: construct a psf (pmPSF) from this collection of model fits
     pmPSFFromModels (test->psf, test->modelFLT, test->mask);
+    
+    // count valid sources
+    int Nkeep = 0;
+    for (int i = 0; i < sources->n; i++) {
+      if (test->mask->data.U8[i]) continue;
+      Nkeep++;
+    }
+    psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates\n", Nkeep, sources->n);
 
     // stage 3: refit with fixed shape parameters
@@ -96,4 +105,5 @@
     }
     psLogMsg ("psphot.psftest", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
+    psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
     DumpModelFits (test->modelPSF, "modelsPSF.dat");
 
@@ -152,6 +162,4 @@
 bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask) {
 
-    int n;
-
     // construct the fit vectors from the collection of objects
     // use the mask to ignore missing fits 
@@ -187,4 +195,13 @@
 	psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
 	// psPolynomial2DDump (psf->params->data[i]);
+
+	// count valid sources
+	int Nkeep = 0;
+	for (int j = 0; j < mask->n; j++) {
+	  if (mask->data.U8[j]) continue;
+	  Nkeep++;
+	}
+	psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
+
     }
     return (true);
Index: /trunk/psphot/src/source_moments.c
===================================================================
--- /trunk/psphot/src/source_moments.c	(revision 4128)
+++ /trunk/psphot/src/source_moments.c	(revision 4129)
@@ -18,6 +18,4 @@
     psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
     keep           = psRegionForImage (keep, image, keep);
-    peaks          = pmPeaksSubset (allpeaks, 100000, keep);
-
     // I need to exclude the saturated objects here, but include them again later...?
     peaks          = pmPeaksSubset (allpeaks, 100000, keep);
Index: /trunk/psphot/src/subtract_galaxies.c
===================================================================
--- /trunk/psphot/src/subtract_galaxies.c	(revision 4128)
+++ /trunk/psphot/src/subtract_galaxies.c	(revision 4129)
@@ -5,6 +5,4 @@
 bool subtract_galaxies (psArray *sources, psMetadata *config) 
 { 
-    bool  status;
-
     for (int i = 0; i < sources->n; i++) {
 	psSource *source = sources->data[i];
Index: /trunk/psphot/src/subtract_psf_source.c
===================================================================
--- /trunk/psphot/src/subtract_psf_source.c	(revision 4128)
+++ /trunk/psphot/src/subtract_psf_source.c	(revision 4129)
@@ -4,8 +4,7 @@
 // need to control application of this with some thresholding
 
-bool subtract_psf_sources (psSource *source)
+bool subtract_psf_source (psSource *source)
 { 
   float sky;
-  psSource *source = sources->data[i];
 
   // non-stellar sources are ignored
@@ -15,5 +14,5 @@
   psModel *model = source->modelPSF;
   if (model == NULL) {
-    psLogMsg ("psphot.subtract_psf_source", "missing model for %f, %f\n", source->moments->x, source->moments->y);
+    psLogMsg ("psphot.subtract_psf_source", 2, "missing model for %f, %f\n", source->moments->x, source->moments->y);
     return (false); 
   }	    
Index: /trunk/psphot/src/test_psf_scatter.c
===================================================================
--- /trunk/psphot/src/test_psf_scatter.c	(revision 4128)
+++ /trunk/psphot/src/test_psf_scatter.c	(revision 4129)
@@ -70,8 +70,9 @@
     snBin ->n = bin;
 
+    // XXX these APIs should get changed
     psPolynomial1D *dsxLine = psPolynomial1DAlloc (2, PS_POLYNOMIAL_ORD);
-    psVectorFitPolynomial1D (dsxLine, NULL, snBin, dsxBin, NULL);
+    psVectorFitPolynomial1D (dsxLine, snBin, dsxBin, NULL);
     psPolynomial1D *dsyLine = psPolynomial1DAlloc (2, PS_POLYNOMIAL_ORD);
-    psVectorFitPolynomial1D (dsyLine, NULL, snBin, dsyBin, NULL);
+    psVectorFitPolynomial1D (dsyLine, snBin, dsyBin, NULL);
     // these should have a slope of 1.0
     psPolynomial1DDump (dsxLine);