Changeset 30883

Timestamp:

Mar 12, 2011, 10:06:53 AM (15 years ago)

Author:

eugene

Message:

limit the number of psf clump regions based on number of available sources; require better sources in the initial assessment of the window for the psf clump (SNLIM*2)

Location:

branches/eam_branches/ipp-20110213/psphot/src

Files:

• psphot.h (modified) (1 diff)
• psphotMergeSources.c (modified) (3 diffs)
• psphotRoughClass.c (modified) (2 diffs)
• psphotSourceStats.c (modified) (6 diffs)
• psphotVisual.c (modified) (1 diff)

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

@@ -253 +253 @@
 bool            psphotVisualShowLogSignificance (psImage *image, float min, float max);
 bool            psphotVisualShowPeaks (pmDetections *detections);
+bool            psphotVisualShowSources (psArray *sources);
 bool            psphotVisualShowFootprints (pmDetections *detections);
 bool            psphotVisualShowMoments (psArray *sources);

branches/eam_branches/ipp-20110213/psphot/src/psphotMergeSources.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
+// Mask to apply for PSF sources : only exclude bad sources -- we will re-test for extendedness
 #define PSF_SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
                          PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
-                         PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply for PSF sources
+                         PM_SOURCE_MODE_CR_LIMIT)
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -347 +348 @@
 
     psLogMsg ("psphot", 3, "%ld PSF sources loaded", detections->peaks->n);
+    psphotVisualShowSources (sources);
+    psphotVisualShowPeaks (detections);
 
     // save detections on the readout->analysis
@@ -354 +357 @@
     }
     psFree (detections);
+
     return true;
 }
 
 // generate the detection structure for the supplied array of sources
+// XXX this function is currently unused
 bool psphotSetSourceParams (pmConfig *config, psArray *sources, pmPSF *psf) {
 

branches/eam_branches/ipp-20110213/psphot/src/psphotRoughClass.c

@@ -63 +63 @@
     }
 
+    int NstarsInClump = psMetadataLookupS32 (&status, readout->analysis, "PSF_CLUMP_NSTARS");
+    // if NstarsInClump is not defined, use the user-selected option:
+    if (!status) {
+        NstarsInClump = 1000;
+    }
+
+    int ScaleForClump = 1;
+    if (NstarsInClump >= 12) ScaleForClump = 2; // 4 cells
+    if (NstarsInClump >= 27) ScaleForClump = 3; // 9 cells
+    if (NstarsInClump >= 48) ScaleForClump = 4; // 16 cells
+    if (NstarsInClump >  75) ScaleForClump = 5; // 25 cells
+
     // we make this measurement on a NxM grid of regions across the readout
     int NX  = psMetadataLookupS32 (&status, recipe, "PSF_CLUMP_NX");  CHECK_STATUS (status, "PSF_CLUMP_NX");
     int NY  = psMetadataLookupS32 (&status, recipe, "PSF_CLUMP_NY");  CHECK_STATUS (status, "PSF_CLUMP_NY");
-    int dX  = readout->image->numCols / NX;
-    int dY  = readout->image->numRows / NY;
+
+    int NXuse, NYuse;
+
+    ScaleForClump = PS_MIN(ScaleForClump, PS_MAX(NX, NY));
+    if (NX > NY) {
+        NXuse = ScaleForClump;
+        NYuse = (int) (ScaleForClump * (NY / NX) + 0.5);
+    } else {
+        NYuse = ScaleForClump;
+        NXuse = (int) (ScaleForClump * (NY / NX) + 0.5);
+    }
+
+    psLogMsg ("psphot", 4, "With %d stars, using %d x %d grid for PSF clump\n", NstarsInClump, NXuse, NYuse);
+
+    int dX  = readout->image->numCols / NXuse;
+    int dY  = readout->image->numRows / NYuse;
 
     int nRegion = 0;
-    for (int ix = 0; ix < NX; ix ++) {
-        for (int iy = 0; iy < NY; iy ++) {
+    for (int ix = 0; ix < NXuse; ix ++) {
+        for (int iy = 0; iy < NYuse; iy ++) {
 
             psRegion *region = psRegionAlloc (ix*dX, (ix + 1)*dX, iy*dY, (iy + 1)*dY);
@@ -181 +207 @@
         return false;
     }
-    psLogMsg ("psphot", 3, "psf clump  X,  Y: %f, %f\n", psfClump.X, psfClump.Y);
-    psLogMsg ("psphot", 3, "psf clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY);
+    psLogMsg ("psphot", 3, "psf clump  X,  Y: %f, %f : DX, DY: %f, %f : nStars %d of %d\n", psfClump.X, psfClump.Y, psfClump.dX, psfClump.dY, psfClump.nStars, psfClump.nTotal);
 
     // get basic parameters, or set defaults

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

@@ -482 +482 @@
 
     float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); psAssert (status, "missing MOMENTS_SN_MIN");
-    float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
     psF32 MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX");
+
+    // when we set the window, we are not attempting to measure spatial variations; we can use a somewhat higher S/N limit
+    // since we are using all sources (true?)
+    float PSF_SN_LIM = 2.0*psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
 
     // XXX move this to a config file?
@@ -489 +492 @@
     float sigma[NSIGMA] = {1.0, 2.0, 3.0, 4.5, 6.0, 9.0, 12.0, 18.0};
     float Sout[NSIGMA];
+    int   Nout[NSIGMA]; // number of stars found in clump : use this to control the number of regions measured by psphotRoughClass
 
     // this sorts by peak->SN
@@ -519 +523 @@
         // 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_AR_MAX);
-        psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
+        psLogMsg ("psphot", 3, "radius %.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]);
 
 #if 0
@@ -540 +544 @@
 
         Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
+        Nout[i] = psfClump.nStars;
     }
 
     // we are looking for sigma for which Sout = 0.65 (or some other value)
+    int Nstars = 0;
     float Sigma = NAN;
     float minS = Sout[0];
@@ -558 +564 @@
         if ((Sout[i] < 0.65) && (Sout[i+1] < 0.65)) continue;
         Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]);
+        Nstars = 0.5*(Nout[i] + Nout[i+1]);
     }
     psAssert (isfinite(Sigma), "did we miss a case?");
@@ -567 +574 @@
     psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma);
     psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
+    psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_CLUMP_NSTARS", PS_META_REPLACE, "number of stars in clump", Nstars);
 
     psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);

branches/eam_branches/ipp-20110213/psphot/src/psphotVisual.c

@@ -350 +350 @@
 }
 
+// requires psphotVisualShowImage
+bool psphotVisualShowSources (psArray *sources) {
+
+    int Noverlay;
+    KiiOverlay *overlay;
+
+    if (!pmVisualTestLevel("psphot.objects.sources", 1)) return true;
+
+    int kapa = psphotKapaChannel (1);
+    if (kapa == -1) return false;
+
+    // note: this uses the Ohana allocation tools:
+    // ALLOCATE (overlay, KiiOverlay, 3*peaks->n + 1);
+    ALLOCATE (overlay, KiiOverlay, sources->n + 2);
+
+    Noverlay = 0;
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+        if (!source) continue;
+
+        pmPeak *peak = source->peak;
+        if (!peak) continue;
+
+        overlay[Noverlay].type = KII_OVERLAY_BOX;
+        overlay[Noverlay].x = peak->xf;
+        overlay[Noverlay].y = peak->yf;
+        overlay[Noverlay].dx = 4.0;
+        overlay[Noverlay].dy = 4.0;
+        overlay[Noverlay].angle = 0.0;
+        overlay[Noverlay].text = NULL;
+        Noverlay ++;
+    }
+
+    KiiLoadOverlay (kapa, overlay, Noverlay, "blue");
+    FREE (overlay);
+
+    pmVisualAskUser(NULL);
+    return true;
+}
+
 // XXX : requires psphotVisualShowImage
 bool psphotVisualShowPeaks (pmDetections *detections) {