Changeset 34415

Timestamp:

Sep 7, 2012, 10:06:48 AM (14 years ago)

Author:

eugene

Message:

when growing the CR mask, copy the values into the existing readout->mask array, otherwise we lose the association to source->maskView; mask a small region for sources that have already been detected (the exclusion radius is curretnly set to 4 pixels, and is not user-specified; deactivate some test code; saturate linear fit fluxes to min valid flux (probably superfluous); calls to psphotFitSourcesLinear in psphotReadout should skip negative fluxes

Location:

branches/eam_branches/ipp-20120905/psphot/src

Files:

• psphotAddNoise.c (modified) (6 diffs)
• psphotBlendFit.c (modified) (5 diffs)
• psphotDeblendSatstars.c (modified) (2 diffs)
• psphotFindDetections.c (modified) (1 diff)
• psphotFitSourcesLinear.c (modified) (1 diff)
• psphotKronIterate.c (modified) (2 diffs)
• psphotReadout.c (modified) (4 diffs)
• psphotSourceSize.c (modified) (1 diff)
• psphotStackReadout.c (modified) (1 diff)

branches/eam_branches/ipp-20120905/psphot/src/psphotAddNoise.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+
+bool psphotMaskSource(pmSource *source, bool add, psImageMaskType maskVal);
 
 bool psphotAddNoise (pmConfig *config, const pmFPAview *view, const char *filerule) {
@@ -30 +32 @@
 }
 
+static int Nmasked = 0;
+
 bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add) {
 
@@ -56 +60 @@
     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     psAssert (maskVal, "missing mask value?");
+
+    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels
+    psAssert (markVal, "missing mask value?");
 
     // increase variance by factor*(object noise):
@@ -94 +101 @@
 
         pmSourceNoiseOp (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, FACTOR, SIZE, add, maskVal, 0, 0);
+
+        psphotMaskSource (source, add, markVal);
     }
     if (add) {
@@ -100 +109 @@
         psLogMsg ("psphot.noise", PS_LOG_WARN, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
     }
+    fprintf (stderr, "masked %d objects\n", Nmasked);
 
     psphotVisualShowImage (readout);
@@ -105 +115 @@
     return true;
 }
+
+bool psphotMaskSource(pmSource *source, bool add, psImageMaskType maskVal) {
+
+    if (!source) return false;
+    if (!source->peak) return false; // XXX how can we have a peak-less source?
+    if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
+    if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
+
+    float Xc = source->peak->xf - source->pixels->col0 - 0.5;
+    float Yc = source->peak->yf - source->pixels->row0 - 0.5;
+
+    psImageMaskType notMaskVal = ~maskVal;
+
+    for (int iy = 0; iy < source->pixels->numRows; iy++) {
+        for (int ix = 0; ix < source->pixels->numCols; ix++) {
+
+            float radius = hypot (ix - Xc, iy - Yc) ;
+
+            if (radius > 4) continue;
+
+            if (add) {
+              source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= maskVal;
+            } else {
+              source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] &= notMaskVal;
+            }
+        }
+    }
+    Nmasked ++;
+
+    return true;
+}
+

branches/eam_branches/ipp-20120905/psphot/src/psphotBlendFit.c

@@ -96 +96 @@
     assert (status && isfinite(skySig) && skySig > 0);
 
-# if (1)
+# if (0)
     // **** test block : generate an ID image where pixels are set based the source models (flux > 0.1 peak && flux > 2.0 skySig)
     psImage *IDimage = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_S32);
@@ -258 +258 @@
 
         int TEST_ON = false;
-# if (1)
+# if (0)
 # define TEST_X 653
 # define TEST_Y 466
         if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
             fprintf (stderr, "test object\n");
-            //  psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
+            psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
             TEST_ON = true;
         }
@@ -379 +379 @@
     if (!source) return false;
     if (!source->peak) return false; // XXX how can we have a peak-less source?
+
+# if (0)
+# define TEST_X 653
+# define TEST_Y 466
+        if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+            fprintf (stderr, "test object\n");
+        }
+# undef TEST_X
+# undef TEST_Y
+# endif
+
     if (!source->moments) return false;
     if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
@@ -473 +484 @@
 
     // what is the radius for a specific peak fraction?
-    for (int i = 0; i < logFmodel->n; i++) {
+    for (int i = 1; i < logFmodel->n; i++) {
         float logF = logFmodel->data.F32[i];
+        if (!isfinite(logF)) continue;
+
         float flux = pow(10.0, logF);
         if (flux > threshold) continue;
-        if (i == 0) continue;
         
         float logF0 = logFmodel->data.F32[i - 1];
@@ -512 +524 @@
     int ID = source->id;
 
+    fprintf (stderr, "Xo,Yo: %f,%f; radius: %f\n", Xo, Yo, radius);
+
     for (int iy = minY; iy < maxY; iy++) {
         for (int ix = minX; ix < maxX; ix++) {

branches/eam_branches/ipp-20120905/psphot/src/psphotDeblendSatstars.c

@@ -121 +121 @@
     int BIG_SIGMA = BIG_RADIUS / 4.0;
 
-    int display = psphotKapaChannel (1);
-    psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 0);
+    // int display = psphotKapaChannel (1);
+    // psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 0);
 
     // examine sources in decreasing SN order
@@ -182 +182 @@
     }
     // show the image after object have been subtracted
-    psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 1);
+    // psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 1);
 
     psFree (SN);

branches/eam_branches/ipp-20120905/psphot/src/psphotFindDetections.c

@@ -45 +45 @@
     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     psAssert (maskVal, "missing mask value?");
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels
+    psAssert (markVal, "missing mark value?");
+
+    maskVal |= markVal;
 
     // Use the new pmFootprints approach?

branches/eam_branches/ipp-20120905/psphot/src/psphotFitSourcesLinear.c

@@ -439 +439 @@
         model->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
         model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
+
+        if (norm->data.F32[i] < MIN_VALID_FLUX) {
+          fprintf (stderr, "fit out of bounds for %f,%f : %f\n", source->peak->xf, source->peak->yf, norm->data.F32[i]);
+          model->params->data.F32[PM_PAR_I0] = MIN_VALID_FLUX;
+        }
 
         // clear the 'mark' pixels so the subtraction covers the full window

branches/eam_branches/ipp-20120905/psphot/src/psphotKronIterate.c

@@ -321 +321 @@
             if (!source->peak) continue; // XXX how can we have a peak-less source?
 
+# if (0)
+# define TEST_X 653
+# define TEST_Y 466
+        if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+            fprintf (stderr, "test object\n");
+        }
+# undef TEST_X
+# undef TEST_Y
+# endif
+
             // check status of this source's moments
             if (!source->moments) continue;
@@ -541 +551 @@
     if (RF <= 0. || RS <= 0 || !isfinite(MrfTry)) {
         // We did not get a good measurement
-        source->moments->Mrf = NAN;
+        source->moments->Mrf = NAN;
         source->moments->KronFlux  = NAN;
         source->moments->KronFluxErr  = NAN;

branches/eam_branches/ipp-20120905/psphot/src/psphotReadout.c

@@ -185 +185 @@
 
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
-    psphotFitSourcesLinear (config, view, filerule, false, false); // pass 1 (detections->allSources)
+    psphotFitSourcesLinear (config, view, filerule, false, true); // pass 1 (detections->allSources)
 
     // measure the radial profiles to the sky
@@ -212 +212 @@
     // linear fit to include all sources (subtract again)
     // NOTE : apply to ALL sources (extended + psf)
-    psphotFitSourcesLinear (config, view, filerule, true, false); // pass 2 (detections->allSources)
+    psphotFitSourcesLinear (config, view, filerule, true, true); // pass 2 (detections->allSources)
 
     // if we only do one pass, skip to extended source analysis
@@ -260 +260 @@
 
         // NOTE: apply to ALL sources
-        psphotFitSourcesLinear (config, view, filerule, true, false); // pass 3 (detections->allSources)
+        psphotFitSourcesLinear (config, view, filerule, true, true); // pass 3 (detections->allSources)
     }
 
@@ -300 +300 @@
 
         // NOTE: apply to ALL sources
-        psphotFitSourcesLinear (config, view, filerule, true, false); // pass 3 (detections->allSources)
+        psphotFitSourcesLinear (config, view, filerule, true, true); // pass 3 (detections->allSources)
     }
 

branches/eam_branches/ipp-20120905/psphot/src/psphotSourceSize.c

@@ -675 +675 @@
             return false;
         }
-        psFree(readout->mask);
-        readout->mask = newMask;
+        // Copy the new mask pixel values to the old mask array.  NOTE: we cannot replace
+        // the mask pointer because the source->maskView objects point to the old data
+        // area
+        psImage *oldMask = readout->mask;
+        readout->mask = psImageCopy (readout->mask, newMask, PS_TYPE_IMAGE_MASK);
+        psAssert (oldMask == readout->mask, "should have copied the values in-situ");
+        psFree(newMask);
     }
 

branches/eam_branches/ipp-20120905/psphot/src/psphotStackReadout.c

@@ -332 +332 @@
     psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
 
-    if (splitLinearFit) {
-        // NOTE: apply to Matched sources. Since the sources that we fit above are subtracted, they will
-        // not be included in this fit.
-        psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources)
-    } else {
-        // Fit all sources together
-        psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 3 (detections->allSources)
-    }
+    // final linear fit. NOTE: if splitLinearFit is true above, this pass will only fit
+    // the unsubtracted (matched) sources (the sources that we fit above are subtracted)
+    psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources)
 
     // measure the radial profiles to the sky (only measures new objects)