Changeset 42439

Timestamp:

Mar 26, 2023, 11:04:35 AM (3 years ago)

Author:

eugene

Message:

psphotFindFeatures now works well to reject sources on trails from Moments and PSF analysis

Location:

branches/eam_branches/ipp-20230313

Files:

• psModules/src/objects/pmSource.c (modified) (2 diffs)
• psphot/src/psphot.h (modified) (1 diff)
• psphot/src/psphotFindFeatures.c (modified) (7 diffs)
• psphot/src/psphotSourceStats.c (modified) (1 diff)

branches/eam_branches/ipp-20230313/psModules/src/objects/pmSource.c

@@ -456 +456 @@
             }
 
+            // skip sources associated with possible detector features 
+            if (source->mode2 & PM_SOURCE_MODE2_ON_LINE) continue;
+
             float Mxx = source->moments->Mxx, Myy = source->moments->Myy; // Second moments
             float ar = Mxx / Myy;       // Radius
@@ -752 +755 @@
                 continue;
             }
+
+            // skip sources associated with possible detector features 
+            if (source->mode2 & PM_SOURCE_MODE2_ON_LINE) continue;
 
             // the rest are probable stellar objects

branches/eam_branches/ipp-20230313/psphot/src/psphot.h

@@ -592 +592 @@
 extern bool psphotINpsphotStack;
 
-bool psphotFindFeatures (psArray *sources);
+bool psphotFindFeatures (pmReadout *readout, psArray *sources);
 
 

branches/eam_branches/ipp-20230313/psphot/src/psphotFindFeatures.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+
+/* Use Hough transform to detect strong linear features (e.g., satellite trails) in the image.
+   Sources which are plausibly associated with the feature are marked in the MODE2 bits.  These 
+   sources are excluded from the analysis of the moments and selection of the PSF stars.
+   Weaknesses: 
+   - the angular samples and the length bins are hard-wired at 1 degree x 10 pixels.
+   - the minimum significance in the Hough image and the minimum number of sources associated
+     with the feature are hard-wired.
+   - the distance of a source from the feature is only calculated in 1D.  The code is not careful
+     about the center of the linear feature compared to the edges of the bin.
+   - the features detected in the Hough transform are used to select sources which are then
+     fitted to a line.  The line is used to select sources again and re-fit.  Only two steps
+     are rigidly defined without a clear iteration end criterion.
+
+ */
 
 # define DCOS(THETA) cos(THETA*PS_RAD_DEG)
@@ -24 +39 @@
 bool psphotFindFeatures_HoughSigma (psphotFindFeatures_HoughInfo *myHough);
 float psphotFindFeatures_LineLength (float radius, float theta, int *edges, int NXpix, int NYpix);
+int psphotFindFeatures_PeaksOnTrail (psArray *sources, float theta, float radius, float Rdel);
+bool psphotFindFeatures_HoughFree (psphotFindFeatures_HoughInfo *myHough);
 
 bool psphotSavePeaks (char *filename, psArray *sources) {
@@ -38 +55 @@
 }
 
+bool psphotSavePoints (char *filename, psVector *x, psVector *y) {
+
+    FILE *ftest = fopen (filename, "w");
+    for (int i = 0; i < x->n; i++) {
+        fprintf (ftest, "%f %f\n", x->data.F32[i], y->data.F32[i]);
+    }
+    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) {
+bool psphotFindFeatures (pmReadout *readout, psArray *sources) {
     
     psphotFindFeatures_HoughInfo myHough;
     
+    int NXpix = readout->image->numCols;
+    int NYpix = readout->image->numRows;
+
     // calculate the Hough transform for the source positions
-    psphotFindFeatures_HoughInit (&myHough, 5000, 5000, 10.0);
-    psphotSaveImage (NULL, myHough.fArea, "hough.fr.fits");
+    psphotFindFeatures_HoughInit (&myHough, NXpix, NYpix, 10.0);
     
     // generate the hough transform for the collection of sources
@@ -83 +112 @@
             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 ++;
-                    }
-                }
-            }
-            
+            psLogMsg ("psphot", PS_LOG_INFO, "trail %d : %f %f (%f : %f)\n", i, theta, radius, hough, nsigma);
+
+            int nPeaksOnThisTrail = psphotFindFeatures_PeaksOnTrail (sources, theta, radius, myHough.Rdel);
+
             if (nPeaksOnThisTrail < HOUGH_MIN_SOURCES) continue;
 
@@ -123 +125 @@
         psFree (lines);
     }
-    fprintf (stderr, "found %d peaks in %d trails\n", nPeaksOnTrails, nTrails);
-    psphotSavePeaks ("test.peaks.dat", sources);
-
+    psphotFindFeatures_HoughFree (&myHough);
+    psLogMsg ("psphot", PS_LOG_INFO, "found %d peaks in %d trails\n", nPeaksOnTrails, nTrails);
+
+    return true;
+}
+
+bool psphotFindFeatures_HoughFree (psphotFindFeatures_HoughInfo *myHough) {
+
+    psFree (myHough->fArea);
+    psFree (myHough->image);
+    psFree (myHough->sigma);
     return true;
 }
@@ -221 +231 @@
     }
     return true;
+}
+
+// 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.
+
+int psphotFindFeatures_PeaksOnTrail (psArray *sources, float theta, float radius, float Rdel) {
+
+    float cosTheta = DCOS(theta);
+    float sinTheta = DSIN(theta);
+
+    // find sources in the line and fit to 1D polynomial
+    psVector *xpos = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *ypos = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
+
+    // if the angle is shallow, use image (x,y), otherwise use image (y,x)
+    bool noFlip = !((theta > 45) && (theta < 135));
+
+    // select only sources close to the line
+    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 (noFlip) {
+            float yLine = (radius - peak->x*cosTheta) / sinTheta;
+            if (fabs(peak->y - yLine) > Rdel) continue;
+            psVectorAppend (xpos, peak->xf);
+            psVectorAppend (ypos, peak->yf);
+        } else {
+            float xLine = (radius - peak->y*sinTheta) / cosTheta;
+            if (fabs(peak->x - xLine) > Rdel) continue;
+            psVectorAppend (xpos, peak->yf);
+            psVectorAppend (ypos, peak->xf);
+        }
+    }
+    if (xpos->n < HOUGH_MIN_SOURCES) goto escape;
+
+    // fit the line
+    if (!psVectorFitPolynomial1D (poly, NULL, 0, ypos, NULL, xpos)) {
+        psError(PS_ERR_UNKNOWN, true, "Failed to find a good trail fit");
+        goto escape;
+    }
+
+    // reset vector lengths
+    xpos->n = 0;
+    ypos->n = 0;
+
+    // select sources close to the line and refit
+    // select only sources close to the line
+    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 (noFlip) {
+            float yLine = psPolynomial1DEval (poly, peak->x);
+            if (fabs(peak->y - yLine) > Rdel) continue;
+            psVectorAppend (xpos, peak->xf);
+            psVectorAppend (ypos, peak->yf);
+        } else {
+            float xLine = psPolynomial1DEval (poly, peak->y);
+            if (fabs(peak->x - xLine) > Rdel) continue;
+            psVectorAppend (xpos, peak->yf);
+            psVectorAppend (ypos, peak->xf);
+        }
+    }
+    if (xpos->n < HOUGH_MIN_SOURCES) goto escape;
+
+    if (!psVectorFitPolynomial1D (poly, NULL, 0, ypos, NULL, xpos)) {
+        psError(PS_ERR_UNKNOWN, true, "Failed to find a good trail fit");
+        goto escape;
+    }
+
+    // mark sources close to the line
+    int nPeaksOnThisTrail = 0;
+    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 (noFlip) {
+            float yLine = psPolynomial1DEval (poly, peak->x);
+            if (fabs(peak->y - yLine) > Rdel) continue;
+        } else {
+            float xLine = psPolynomial1DEval (poly, peak->y);
+            if (fabs(peak->x - xLine) > Rdel) continue;
+        }
+        source->mode2 |= PM_SOURCE_MODE2_ON_LINE;
+        nPeaksOnThisTrail ++;
+    }
+    if (nPeaksOnThisTrail < HOUGH_MIN_SOURCES) goto escape;
+
+    psFree (xpos);
+    psFree (ypos);
+    psFree (poly);
+    return nPeaksOnThisTrail;
+
+escape:
+    psFree (xpos);
+    psFree (ypos);
+    psFree (poly);
+    return false;
 }
 
@@ -320 +435 @@
     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 ++;
-                    }
-                }
-            }
-        }
-*/

branches/eam_branches/ipp-20230313/psphot/src/psphotSourceStats.c

@@ -144 +144 @@
     }
 
-    // XXX I want to identify sources which are on lines.  these can be excluded from
-    // the Moments and PSF analysis
-
-    psphotFindFeatures (sources);
-
     if (setWindow) {
+        // identify sources on lines, drop from MomentsWindow analysis
+        psphotFindFeatures (readout, sources);
+
         if (!psphotSetMomentsWindow(recipe, readout->analysis, sources, maskVal)) {
             psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");