Index: /branches/eam_branches/ipp-20230313/psphot/src/Makefile.am
===================================================================
--- /branches/eam_branches/ipp-20230313/psphot/src/Makefile.am	(revision 42434)
+++ /branches/eam_branches/ipp-20230313/psphot/src/Makefile.am	(revision 42435)
@@ -247,4 +247,5 @@
 	psphotSetNFrames.c	       \
 	psphotSourceMemory.c	       \
+	psphotFindFeatures.c	       \
 	psphotChipParams.c		\
 	psphotGalaxyParams.c
Index: /branches/eam_branches/ipp-20230313/psphot/src/psphot.h
===================================================================
--- /branches/eam_branches/ipp-20230313/psphot/src/psphot.h	(revision 42434)
+++ /branches/eam_branches/ipp-20230313/psphot/src/psphot.h	(revision 42435)
@@ -592,4 +592,6 @@
 extern bool psphotINpsphotStack;
 
+bool psphotFindFeatures (psArray *sources);
+
 
 #endif
Index: /branches/eam_branches/ipp-20230313/psphot/src/psphotFindFeatures.c
===================================================================
--- /branches/eam_branches/ipp-20230313/psphot/src/psphotFindFeatures.c	(revision 42435)
+++ /branches/eam_branches/ipp-20230313/psphot/src/psphotFindFeatures.c	(revision 42435)
@@ -0,0 +1,374 @@
+# include "psphotInternal.h"
+
+# define DCOS(THETA) cos(THETA*PS_RAD_DEG)
+# define DSIN(THETA) sin(THETA*PS_RAD_DEG)
+# define HOUGH_MIN_SOURCES 10
+# define HOUGH_MIN_SIGMA 5.0
+
+typedef struct {
+    psImage *fArea;
+    psImage *image;
+    psImage *sigma;
+    int NXpix;
+    int NYpix;
+    float Rdel;
+    float Rmin;
+    float Rmax;
+    int Nt;
+    int Nr;
+    int nKeep;
+} psphotFindFeatures_HoughInfo;
+
+bool psphotFindFeatures_HoughInit (psphotFindFeatures_HoughInfo *myHough, int NXpix, int NYpix, float Rdel);
+bool psphotFindFeatures_HoughImage (psphotFindFeatures_HoughInfo *myHough, psArray *sources);
+bool psphotFindFeatures_HoughSigma (psphotFindFeatures_HoughInfo *myHough);
+float psphotFindFeatures_LineLength (float radius, float theta, int *edges, int NXpix, int NYpix);
+
+bool psphotSavePeaks (char *filename, psArray *sources) {
+
+    FILE *ftest = fopen (filename, "w");
+    for (int i = 0; i < sources->n; i++) {
+	pmSource *source = sources->data[i];
+	pmPeak *peak = source->peak;
+	bool onTrail = source->mode2 & PM_SOURCE_MODE2_ON_LINE;
+	fprintf (ftest, "%d %d %f %d\n", peak->x, peak->y, peak->detValue, onTrail);
+    }
+    fclose (ftest);
+    return true;
+}
+
+// examine the x,y distribution of the sources and attempt to find sources which are due to
+// specific finds of features.  For example, a bunch of detections in a line could be a
+// satellite and should not be used for PSF modeling
+
+bool psphotFindFeatures (psArray *sources) {
+    
+    psphotFindFeatures_HoughInfo myHough;
+    
+    // calculate the Hough transform for the source positions
+    psphotFindFeatures_HoughInit (&myHough, 5000, 5000, 10.0);
+    psphotSaveImage (NULL, myHough.fArea, "hough.fr.fits");
+    
+    // generate the hough transform for the collection of sources
+    // and find only a single trail at a time.  remove that trail
+    // and try again.
+
+    int nTrails = 0;
+    int nPeaksOnTrails = 0;
+    bool foundTrail = true; // true to ensure we try a first pass
+    while (foundTrail) {
+	psphotFindFeatures_HoughImage (&myHough, sources);
+	psphotFindFeatures_HoughSigma (&myHough);
+	
+	// select the peaks, sort by the peak amplitude
+	psArray *lines = pmPeaksInImage (myHough.sigma, HOUGH_MIN_SIGMA);
+	psArraySort (lines, pmPeaksSortByRawFluxDescend);
+	
+	// if we do not find a trail in this pass, we stop
+	foundTrail = false;
+
+	psImage *image = myHough.image;
+	psImage *sigma = myHough.sigma;
+
+	// choose the strongest line (> HOUGH_MIN_SIGMA) with at least HOUGH_MIN_SOURCES
+	// and mark sources associated with this line
+	for (int i = 0; (i < lines->n) && !foundTrail; i++) {
+	    pmPeak *line = lines->data[i];
+
+	    float hough = image->data.F32[line->y][line->x];
+	    if (hough < HOUGH_MIN_SOURCES) continue;
+
+	    float nsigma = sigma->data.F32[line->y][line->x];
+
+	    float theta  = line->x;
+	    float radius = line->y*myHough.Rdel + myHough.Rmin;
+	    fprintf (stderr, "trail %d : %f %f (%f : %f)\n", i, theta, radius, hough, nsigma);
+
+	    float cosTheta = DCOS(theta);
+	    float sinTheta = DSIN(theta);
+
+	    int nPeaksOnThisTrail = 0;
+
+	    // mark all sources within myHough->Rdel pix of this line
+	    // I should probably calculate the distance from the line in 2D perpendicular direction
+	    // NOTE: I am marking the sources on this trail, but if it does not have enough
+	    // matches, it is not counted as a real peak.
+	    for (int j = 0; j < sources->n; j++) {
+		pmSource *source = sources->data[j];
+		if (source->mode2 & PM_SOURCE_MODE2_ON_LINE) continue; // ignore sources already found
+
+		pmPeak *peak = source->peak;
+		if ((theta > 45) && (theta < 135)) {
+		    float yLine = (radius - peak->x*cosTheta) / sinTheta;
+		    if (fabs(peak->y - yLine) < myHough.Rdel) {
+			source->mode2 |= PM_SOURCE_MODE2_ON_LINE;
+			nPeaksOnThisTrail ++;
+		    }
+		} else {
+		    float xLine = (radius - peak->y*sinTheta) / cosTheta;
+		    if (fabs(peak->x - xLine) < myHough.Rdel) {
+			source->mode2 |= PM_SOURCE_MODE2_ON_LINE;
+			nPeaksOnThisTrail ++;
+		    }
+		}
+	    }
+	    
+	    if (nPeaksOnThisTrail < HOUGH_MIN_SOURCES) continue;
+
+	    // this is only a valid trail if it has enough sources
+	    nTrails ++;
+	    foundTrail = true;
+	    nPeaksOnTrails += nPeaksOnThisTrail;
+	}
+	psFree (lines);
+    }
+    fprintf (stderr, "found %d peaks in %d trails\n", nPeaksOnTrails, nTrails);
+    psphotSavePeaks ("test.peaks.dat", sources);
+
+    return true;
+}
+
+// define the basic parameters for all Hough transform images 
+// the transform image is not populated here, but the fractional area image is
+bool psphotFindFeatures_HoughInit (psphotFindFeatures_HoughInfo *myHough, int NXpix, int NYpix, float Rdel) {
+
+    // diagonal length (absolute max value of radius)
+    int NRpix = sqrt(PS_SQR(NXpix) + PS_SQR(NYpix));
+
+    // radial bins
+    myHough->Rmin = -1*NRpix;
+    myHough->Rmax = +1*NRpix;
+    myHough->Rdel = Rdel;
+
+    myHough->Nt = 180; // angle range from 0 to 180 in steps of 1 degree
+    myHough->Nr = (myHough->Rmax - myHough->Rmin) / myHough->Rdel;
+
+    myHough->NXpix = NXpix;
+    myHough->NYpix = NYpix;
+
+    myHough->fArea = psImageAlloc (myHough->Nt, myHough->Nr, PS_TYPE_F32);
+    myHough->image = psImageAlloc (myHough->Nt, myHough->Nr, PS_TYPE_F32);
+    myHough->sigma = psImageAlloc (myHough->Nt, myHough->Nr, PS_TYPE_F32);
+
+    float ImageArea = myHough->NXpix * myHough->NYpix;
+    for (int theta = 0; theta < myHough->Nt; theta ++) {
+	for (int rbin = 0; rbin < myHough->Nr; rbin ++) {
+
+	    // for the given value of Radius, Theta, what is the line length from edge to edge
+	    int edges[4];
+	    float radius = rbin * myHough->Rdel + myHough->Rmin;
+	    float LineLength = psphotFindFeatures_LineLength (radius, theta, edges, myHough->NXpix, myHough->NYpix);
+	    float LineArea = LineLength * myHough->Rdel;
+	    myHough->fArea->data.F32[rbin][theta] = LineArea / ImageArea;
+	}
+    }
+    return true;
+}
+
+bool psphotFindFeatures_HoughImage (psphotFindFeatures_HoughInfo *myHough, psArray *sources) {
+
+    if (!sources) return false;
+
+    psImage *image = myHough->image;
+    psImageInit (image, 0.0);
+
+    int nKeep = 0;
+
+    // generate the Hough transform image
+    for (int i = 0; i < sources->n; i++) {
+	pmSource *source = sources->data[i];
+	pmPeak *peak = source->peak;
+
+	// skip bad peaks and record the number actually used
+	if (source->mode2 & PM_SOURCE_MODE2_ON_LINE) continue;
+	nKeep ++;
+
+	for (int theta = 0; theta < myHough->Nt; theta ++) {
+	    
+	    float r = peak->x * DCOS(theta) + peak->y * DSIN(theta);
+	    int rbin = (r - myHough->Rmin) / myHough->Rdel;
+
+	    if (theta < 0) continue;
+	    if (rbin  < 0) continue;
+
+	    if (theta >= myHough->Nt) continue;
+	    if (rbin  >= myHough->Nr) continue;
+
+	    image->data.F32[rbin][theta] ++;
+	}
+    }
+    myHough->nKeep = nKeep;
+    return true;
+}
+
+bool psphotFindFeatures_HoughSigma (psphotFindFeatures_HoughInfo *myHough) {
+
+    // Given the fractional area image and nKeep sources
+    // calculate the expected mean number of sources for each line at
+    // radius R and angle theta.  use that to determine significance
+
+    psImage *sigma = myHough->sigma;
+    psImage *image = myHough->image;
+    psImage *fArea = myHough->fArea;
+
+    for (int theta = 0; theta < myHough->Nt; theta ++) {
+	for (int rbin = 0; rbin < myHough->Nr; rbin ++) {
+	    float mean = myHough->nKeep * fArea->data.F32[rbin][theta];
+	    float meanMin = PS_MAX(mean, 1);
+	    float nsigma = (image->data.F32[rbin][theta] - meanMin) / sqrt(meanMin);
+	    sigma->data.F32[rbin][theta] = nsigma;
+	}
+    }
+    return true;
+}
+
+enum {EDGE_LEFT, EDGE_BOTTOM, EDGE_RIGHT, EDGE_TOP};
+
+void psphotFindFeatures_WhichEdges (float radius, float sinTheta, float cosTheta, int *edges, int NXpix, int NYpix) {
+
+    edges[EDGE_TOP]    = false;
+    edges[EDGE_LEFT]   = false;
+    edges[EDGE_RIGHT]  = false;
+    edges[EDGE_BOTTOM] = false;
+    
+    // r = 0 goes through (0,0) -> left edge
+    if ((radius >= 0) && (radius < NYpix*sinTheta)) {
+	edges[EDGE_LEFT] = true;
+    }
+  
+    // if (radius = NXpix*cosTheta) : goes through (Nx,0) -> bottom edge
+    if (cosTheta > 0.0) {
+	if ((radius > 0) && (radius <= NXpix*cosTheta)) {
+	    edges[EDGE_BOTTOM] = true;
+	}
+    } else {
+	if ((radius < 0) && (radius >= NXpix*cosTheta)) {
+	    edges[EDGE_BOTTOM] = true;
+	}
+    }
+
+    // if (radius = NXpix*cosTheta + NYpix*sinTheta) : goes through (Nx,Ny) -> right edge
+    if ((radius > NXpix*cosTheta) && (radius <= NXpix*cosTheta + NYpix*sinTheta)) {
+	edges[EDGE_RIGHT] = true;
+    }
+    
+    // if (radius = NYpix*sinTheta) : goes through (0,Ny) -> top edge
+    if (cosTheta > 0.0) {
+	if ((radius >= NYpix*sinTheta) && (radius < NXpix*cosTheta + NYpix*sinTheta)) {
+	    edges[EDGE_TOP] = true;
+	}
+    } else {
+	if ((radius <= NYpix*sinTheta) && (radius > NXpix*cosTheta + NYpix*sinTheta))
+	    edges[EDGE_TOP] = true;
+    }
+}
+
+float psphotFindFeatures_LineLength (float radius, float theta, int *edges, int NXpix, int NYpix) {
+
+    float cosTheta = DCOS(theta);
+    float sinTheta = DSIN(theta);
+
+    psphotFindFeatures_WhichEdges (radius, sinTheta, cosTheta, edges, NXpix, NYpix);
+
+    if (fabs(theta -   0.0) < 0.01) return ((radius >= 0) && (radius < +NXpix)) ? NYpix : 0.0;
+    if (fabs(theta - 180.0) < 0.01) return ((radius <= 0) && (radius > -NXpix)) ? NYpix : 0.0;
+    if (fabs(theta -  90.0) < 0.01) return ((radius >= 0) && (radius < +NYpix)) ? NXpix : 0.0;
+	
+    if (edges[EDGE_LEFT] && edges[EDGE_RIGHT]) {
+	float yLT = (radius -     0*cosTheta) / sinTheta;
+	float yRT = (radius - NXpix*cosTheta) / sinTheta;
+	float dY = yLT - yRT;
+	float lineLength = hypot(NXpix, dY);
+	return lineLength;
+    }
+    if (edges[EDGE_LEFT] && edges[EDGE_BOTTOM]) {
+	float yLT = (radius -     0*cosTheta) / sinTheta;
+	float xBT = (radius -     0*sinTheta) / cosTheta;
+	float lineLength = hypot(xBT, yLT);
+	return lineLength;
+    }
+    if (edges[EDGE_LEFT] && edges[EDGE_TOP]) {
+	float yLT = (radius -     0*cosTheta) / sinTheta;
+	float xTP = (radius - NYpix*sinTheta) / cosTheta;
+	float dY = NYpix - yLT;
+	float lineLength = hypot(xTP, dY);
+	return lineLength;
+    }
+
+    if (edges[EDGE_BOTTOM] && edges[EDGE_TOP]) {
+	float xBT = (radius -     0*sinTheta) / cosTheta;
+	float xTP = (radius - NYpix*sinTheta) / cosTheta;
+	float dX = xTP - xBT;
+	float lineLength = hypot(dX, NYpix);
+	return lineLength;
+    }
+    if (edges[EDGE_BOTTOM] && edges[EDGE_RIGHT]) {
+	float xBT = (radius -     0*sinTheta) / cosTheta;
+	float yRT = (radius - NXpix*cosTheta) / sinTheta;
+	float dX = NXpix - xBT;
+	float lineLength = hypot(dX, yRT);
+	return lineLength;
+    }
+    if (edges[EDGE_TOP] && edges[EDGE_RIGHT]) {
+	float xTP = (radius - NYpix*sinTheta) / cosTheta;
+	float yRT = (radius - NXpix*cosTheta) / sinTheta;
+	float dX = NXpix - xTP;
+	float dY = NYpix - yRT;
+	float lineLength = hypot(dX, dY);
+	return lineLength;
+    }
+    return 0.0;
+}
+
+
+/*
+  x = (radius - y*sinTheta) / cosTheta;
+  y = (radius - x*cosTheta) / sinTheta;
+*/
+
+/*
+    // find peaks in this image
+    // what is a statistically significant peak?
+    psStats *stats = psStatsAlloc (PS_STAT_MAX | PS_STAT_SAMPLE_STDEV);
+    if (!psImageStats (stats, HoughImage, NULL, 0)) {
+	psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
+	psFree(stats);
+	psFree(HoughImage);
+	return false;
+    }
+*/
+
+/*
+	int nTrails = 0;
+	for (int iLine = 0; iLine < lines->n; iLine++) {
+	    pmPeak *line = lines->data[iLine];
+	    float theta  = line->x;
+	    float radius = line->y*Rdel + Rmin;
+
+	    float hough = HoughImage->data.F32[line->y][line->x];
+	    float sigma = HoughSigma->data.F32[line->y][line->x];
+	    fprintf (stderr, "trail %d : %f %f (%f : %f)\n", iLine, theta, radius, hough, sigma);
+
+	    float cosTheta = DCOS(theta);
+	    float sinTheta = DSIN(theta);
+
+	    // XXX I can probably calculate the distance from the line in a more exact fashion...
+	    for (int i = 0; i < sources->n; i++) {
+		pmSource *source = sources->data[i];
+		pmPeak *peak = source->peak;
+		if ((theta > 45) && (theta < 135)) {
+		    float yLine = (radius - peak->x*cosTheta) / sinTheta;
+		    if (fabs(peak->y - yLine) < 10.0) {
+			source->mode2 |= PM_SOURCE_MODE2_ON_LINE;
+			nTrails ++;
+		    }
+		} else {
+		    float xLine = (radius - peak->y*sinTheta) / cosTheta;
+		    if (fabs(peak->x - xLine) < 10.0) {
+			source->mode2 |= PM_SOURCE_MODE2_ON_LINE;
+			nTrails ++;
+		    }
+		}
+	    }
+	}
+*/
Index: /branches/eam_branches/ipp-20230313/psphot/src/psphotSourceStats.c
===================================================================
--- /branches/eam_branches/ipp-20230313/psphot/src/psphotSourceStats.c	(revision 42434)
+++ /branches/eam_branches/ipp-20230313/psphot/src/psphotSourceStats.c	(revision 42435)
@@ -120,4 +120,5 @@
         source->moments = pmMomentsAlloc();
 
+	// XXX can this ever be set at this point? (we just allocated the source and moments)
         if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) {
 	    fprintf (stderr, "moment failure\n");
@@ -142,4 +143,9 @@
         return true;
     }
+
+    // XXX I want to identify sources which are on lines.  these can be excluded from
+    // the Moments and PSF analysis
+
+    psphotFindFeatures (sources);
 
     if (setWindow) {
@@ -539,11 +545,23 @@
     sources = psArraySort (sources, pmSourceSortByFlux);
 
+    // generate an array of the sources which we want to use for this analysis
+    // XXX move max source number to config
+    psArray *sourceSubset = psArrayAllocEmpty (100);
+    for (int i = 0; (i < sources->n) && (i < 400); i++) {
+
+	pmSource *source = sources->data[i];
+
+	// skip faint sources for moments measurement
+	if (sqrt(source->peak->detValue) < MIN_SN) continue;
+
+	// skip sources on lines
+	if (source->mode2 & PM_SOURCE_MODE2_ON_LINE) continue;
+	psArrayAdd (sourceSubset, 100, source);
+    }
+
     // loop over radii:
     for (int i = 0; i < nsigma; i++) {
-
-        // XXX move max source number to config
-        for (int j = 0; (j < sources->n) && (j < 400); j++) {
-
-            pmSource *source = sources->data[j];
+        for (int j = 0; j < sourceSubset->n; j++) {
+            pmSource *source = sourceSubset->data[j];
             psAssert (source->moments, "force moments to exist");
             source->moments->nPixels = 0;
@@ -565,6 +583,6 @@
 	  sprintf (name, "moments.v%d.dat", i);
 	  FILE *fout = fopen (name, "w");
-	  for (int j = 0; j < sources->n; j++) {
-            pmSource *source = sources->data[j];
+	  for (int j = 0; j < sourceSubset->n; j++) {
+            pmSource *source = sourceSubset->data[j];
             psAssert (source->moments, "force moments to exist");
             source->moments->nPixels = 0;
@@ -581,8 +599,8 @@
 
         // determine the PSF parameters from the source moment values
-        pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_SX_MIN, MOMENTS_SY_MIN, MOMENTS_AR_MAX);
+        pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sourceSubset, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_SX_MIN, MOMENTS_SY_MIN, MOMENTS_AR_MAX);
         psLogMsg ("psphot", 3, "sigma guess (pix) %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
 
-	Rmin[i] = pmSourceMinKronRadius(sources, PSF_SN_LIM);
+	Rmin[i] = pmSourceMinKronRadius(sourceSubset, PSF_SN_LIM);
 
 #if 0
@@ -600,5 +618,5 @@
         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
 	if (pmVisualTestLevel("psphot.moments.full", 2)) {
-	    psphotVisualPlotMoments (recipe, analysis, sources);
+	    psphotVisualPlotMoments (recipe, analysis, sourceSubset);
 	}
 #endif
@@ -608,4 +626,6 @@
         Sout[i] = (Nout[i] == 0) ? NAN : sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
     }
+
+    psFree (sourceSubset);
 
     // we are looking for sigma for which Sout = 0.65 (or some other value)