Index: /branches/eam_branches/ipp-20230313/psModules/src/objects/pmSource.c =================================================================== --- /branches/eam_branches/ipp-20230313/psModules/src/objects/pmSource.c (revision 42438) +++ /branches/eam_branches/ipp-20230313/psModules/src/objects/pmSource.c (revision 42439) @@ -456,4 +456,7 @@ } + // 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,4 +755,7 @@ continue; } + + // skip sources associated with possible detector features + if (source->mode2 & PM_SOURCE_MODE2_ON_LINE) continue; // the rest are probable stellar objects Index: /branches/eam_branches/ipp-20230313/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20230313/psphot/src/psphot.h (revision 42438) +++ /branches/eam_branches/ipp-20230313/psphot/src/psphot.h (revision 42439) @@ -592,5 +592,5 @@ extern bool psphotINpsphotStack; -bool psphotFindFeatures (psArray *sources); +bool psphotFindFeatures (pmReadout *readout, psArray *sources); Index: /branches/eam_branches/ipp-20230313/psphot/src/psphotFindFeatures.c =================================================================== --- /branches/eam_branches/ipp-20230313/psphot/src/psphotFindFeatures.c (revision 42438) +++ /branches/eam_branches/ipp-20230313/psphot/src/psphotFindFeatures.c (revision 42439) @@ -1,3 +1,18 @@ # 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,4 +39,6 @@ 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,15 +55,27 @@ } +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,35 +112,8 @@ 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,7 +125,15 @@ 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,4 +231,109 @@ } 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,55 +435,2 @@ 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 42438) +++ /branches/eam_branches/ipp-20230313/psphot/src/psphotSourceStats.c (revision 42439) @@ -144,10 +144,8 @@ } - // 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!");