Changeset 34418 for trunk

Timestamp:

Sep 7, 2012, 10:51:18 AM (14 years ago)

Author:

eugene

Message:

merge changes from eam_branch/ipp-20120905: add satstar profile to psphotStackReadout; skip negative fluxes in standard psphotReadout; replace mask pixel values when growing the CR mask; mask cores of 1st pass sources to avoid excessively close detections; add SUBTRACT_SATSTAR_PROFILE to psphot recipe

Location:

trunk

Files:

• . (modified) (1 prop)
• Ohana (modified) (1 prop)
• Ohana/src/opihi (modified) (1 prop)
• Ohana/src/opihi/cmd.astro (modified) (1 prop)
• ippMonitor (modified) (1 prop)
• ippScripts/scripts/destreak_restore_camera.pl (modified) (1 prop)
• ippScripts/scripts/ipp_apply_burntool_single.pl (modified) (1 prop)
• ippScripts/scripts/magic_destreak.pl (modified) (1 prop)
• ippToPsps (modified) (1 prop)
• ippTools/share/camtool_find_pendingimfile.sql (modified) (1 prop)
• ippTools/share/chiptool_setimfiletoupdate.sql (modified) (1 prop)
• ippTools/share/pxadmin_create_tables.sql (modified) (1 prop)
• ippTools/share/warptool_towarped.sql (modified) (1 prop)
• ippTools/src (modified) (1 prop)
• ippTools/src/magictool.c (modified) (1 prop)
• ippconfig (modified) (1 prop)
• ippconfig/recipes/psphot.config (modified) (2 diffs)
• ippconfig/recipes/reductionClasses.mdc (modified) (1 prop)
• ppImage/src (modified) (1 prop)
• psModules (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (1 prop)
• psphot/src/psphot.h (modified) (1 diff)
• psphot/src/psphotAddNoise.c (modified) (6 diffs)
• psphot/src/psphotBlendFit.c (modified) (7 diffs)
• psphot/src/psphotDeblendSatstars.c (modified) (7 diffs)
• psphot/src/psphotFindDetections.c (modified) (1 diff)
• psphot/src/psphotFitSourcesLinear.c (modified) (1 diff)
• psphot/src/psphotKronIterate.c (modified) (2 diffs)
• psphot/src/psphotReadout.c (modified) (4 diffs)
• psphot/src/psphotSourceSize.c (modified) (1 diff)
• psphot/src/psphotStackImageLoop.c (modified) (1 prop)
• psphot/src/psphotStackReadout.c (modified) (3 diffs)
• pstamp/scripts (modified) (1 prop)
• psvideophot (modified) (1 prop)

trunk/ippconfig/recipes/psphot.config

@@ -75 +75 @@
 DEBLEND_SKY_NSIGMA                  F32   10.0
 DEBLEND_MIN_SATURATION              F32   40000.0
+
+# subtract radial profiles for extremely bright stars?
+SUBTRACT_SATSTAR_PROFILE            BOOL  F
 
 # parameters to control the selection of the peak in the Sx,Sy plane
@@ -414 +417 @@
   PETROSIAN_FOR_STARS                 BOOL  FALSE
 
+  # subtract radial profiles for extremely bright stars?
+  SUBTRACT_SATSTAR_PROFILE            BOOL  T
+  
   CONSTANT_PHOTOMETRIC_WEIGHTS        BOOL  TRUE 
   POISSON.ERRORS.PHOT.LMM             BOOL  TRUE

trunk/psphot/src/psphot.h

@@ -81 +81 @@
 
 bool            psphotDeblendSatstars (pmConfig *config, const pmFPAview *view, const char *filerule);
-bool            psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
+bool            psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
 
 bool            psphotBasicDeblend (pmConfig *config, const pmFPAview *view, const char *filerule);

trunk/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;
+}
+

trunk/psphot/src/psphotBlendFit.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 bool psphotBlendFitSetSource(pmSource *source, psImage *IDimage, float sigSigma);
+bool psphotBlendSetSourceSatstar (psImage *image, pmSource *source);
+float InterpolateValues (float X0, float Y0, float X1, float Y1, float X);
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -256 +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;
         }
@@ -377 +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;
@@ -385 +398 @@
     if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) {
         // find the radius at which we hit something like 0.1*SATURATION
-        return false;
+        psphotBlendSetSourceSatstar (IDimage, source);
+        return true;
     }
 
@@ -431 +445 @@
     float rMxx = 0.5 / PS_SQR(shape.sx);
     float rMyy = 0.5 / PS_SQR(shape.sy);
-    float Sxy = -1. * shape.sxy;    // factor of -1 is included to match the previous window function
-                                    // implementation. XXX: Is this correct?
+    float Sxy = shape.sxy;    // factor of -1 is included to match the previous window function
+    // implementation. XXX: Is this correct?
 
     int ID = source->id;
@@ -443 +457 @@
 
             float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY;
+            if (z > 2.311) continue;
 
             float f = Io*exp(-z);
@@ -453 +468 @@
     return true;
 }
+
+bool psphotBlendSetSourceSatstar (psImage *IDimage, pmSource *source) {
+
+    float Xo = source->satstar->Xo;
+    float Yo = source->satstar->Yo;
+    psVector *logRmodel = source->satstar->logRmodel;
+    psVector *logFmodel = source->satstar->logFmodel;
+
+    float lPeak = logFmodel->data.F32[0];
+    float fPeak = pow(10.0, lPeak);
+    float threshold = 0.1*fPeak;
+    float logThresh = log10(threshold);
+
+    float radius = NAN;
+
+    // what is the radius for a specific peak fraction?
+    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;
+        
+        float logF0 = logFmodel->data.F32[i - 1];
+        float logR0 = logRmodel->data.F32[i - 1];
+
+        float logF1 = logFmodel->data.F32[i];
+        float logR1 = logRmodel->data.F32[i];
+
+        float logR = InterpolateValues (logF0, logR0, logF1, logR1, logThresh);
+        radius = pow(10.0, logR);
+        break;
+    }
+
+    if (!isfinite(radius)) {
+        for (int i = logFmodel->n - 1; i >= 0; i--) {
+            if (!isfinite(logFmodel->data.F32[i])) continue;
+            float logR = logRmodel->data.F32[i];
+            radius = pow(10.0, logR);
+            break;
+        }
+    }
+
+    if (!isfinite(radius)) return false;
+
+    int Nx = IDimage->numCols;
+    int Ny = IDimage->numRows;
+
+    // region to mask
+    int minX = PS_MIN(PS_MAX(Xo - radius, 0), Nx - 1);
+    int maxX = PS_MIN(PS_MAX(Xo + radius, 0), Nx - 1);
+    int minY = PS_MIN(PS_MAX(Yo - radius, 0), Ny - 1);
+    int maxY = PS_MIN(PS_MAX(Yo + radius, 0), Ny - 1);
+
+    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++) {
+
+            float dX = (ix - Xo);
+            float dY = (iy - Yo);
+            float R = hypot (dX, dY) ;
+            if (R > radius) continue;
+
+            IDimage->data.S32[iy][ix] = ID;
+        }
+    }
+    return true;
+}
+

trunk/psphot/src/psphotDeblendSatstars.c

@@ -16 +16 @@
 bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view, const char *filerule)
 {
+    bool status = false;
+
     int num = psphotFileruleCount(config, filerule);
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    // perform full extended source non-linear fits?
+    if (!psMetadataLookupBool (&status, recipe, "SUBTRACT_SATSTAR_PROFILE")) {
+        psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source fits\n");
+        return true;
+    }
 
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-        if (!psphotDeblendSatstarsReadout (config, view, filerule, i)) {
+        if (!psphotDeblendSatstarsReadout (config, view, filerule, i, recipe)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for %s entry %d", filerule, i);
             return false;
@@ -40 +52 @@
     4) subtract the radial profile
     
+    There is also support code for calculation of magnitudes and add/subtract the profile (a la pmSourceOp)
+
     TBD:
 
-    * function to replace the radial profile for a source
     * recenter the profile (based on cross-correlation / convolution with 1D profile)
-
-    * raise a bit somewhere for these super saturated stars (if they were so modeled)
-    * consider the subtraction bit : when to raise it?
-
  **/
 
-bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex) {
+bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex, psMetadata *recipe) {
 
     int N;
@@ -56 +65 @@
     bool status;
 
-    // XXX disable this function for now
-    return true;
-
     psTimerStart ("psphot.deblend.sat");
 
@@ -78 +84 @@
         return true;
     }
-
-    // select the appropriate recipe information
-    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
-    psAssert (recipe, "missing recipe?");
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
@@ -124 +126 @@
     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
@@ -185 +187 @@
     }
     // 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);
@@ -451 +453 @@
     // XXX do something sensible here if the profile is crap
 
+    // replace an existing profile
+    psFree (source->satstar);
     source->satstar = pmSourceSatstarAlloc();
     source->satstar->Xo = Xo;

trunk/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?

trunk/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

trunk/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;

trunk/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)
     }
 

trunk/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);
     }
 

trunk/psphot/src/psphotStackReadout.c

@@ -169 +169 @@
     }
 
+    // If DET and SRC are different images, subtract radial profiles for the convolved
+    // image first.  The profiles found on the convolved image will be replaced by those
+    // found for the unconvolved image.  Downstream, in psphotFindDetections, we will
+    // replace the the profiles (and re-subtract them) for the detection image, so we want
+    // to keep those versions of the profiles on the sources.
+    if (strcmp(STACK_SRC, STACK_DET)) {
+      // find and subtract radial profile models for saturated stars (XXX change name eventually)
+      if (!psphotDeblendSatstars (config, view, STACK_SRC)) {
+        psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+      }
+    }
+    // find and subtract radial profile models for saturated stars (XXX change name eventually)
+    if (!psphotDeblendSatstars (config, view, STACK_DET)) {
+        psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
+        return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+
     // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
     // only run this on detections from the input images, not chisq image
@@ -206 +224 @@
     // window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
     // but iterates to an appropriately larger size
-        logMemStats("before.kron.1");
+    logMemStats("before.kron.1");
     psphotKronIterate(config, view, STACK_SRC, 1);
-        logMemStats("after.kron.1");
-
+    logMemStats("after.kron.1");
+        
     // identify CRs and extended sources
     psphotSourceSize (config, view, STACK_SRC, true);
@@ -332 +350 @@
     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)