Changeset 4116

Timestamp:

Jun 5, 2005, 8:41:28 PM (21 years ago)

Author:

eugene

Message:

using source masks to select PSF stars

Location:

trunk/psphot/src

Files:

• apply_psf_model.c (modified) (3 diffs)
• mark_psf_source.c (modified) (2 diffs)
• psPolynomials.c (modified) (1 diff)
• pspsf.c (modified) (1 diff)
• subtract_psf_source.c (modified) (1 diff)

trunk/psphot/src/apply_psf_model.c

@@ -1 +1 @@
 # include "psphot.h"
 
-// fit psf model to all objects not already fitted (modelPSF == NULL)
+// fit psf model to all objects 
+// PSFSTAR objects will be refitted
 // run this function to a specific flux limit?
+
 
 bool apply_psf_model (psImage *image, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) 
@@ -35 +37 @@
         psSource *source = sources->data[i];
 
-        // skip the existing PSF stars 
-        // XXX drop this -- refit?
-        if (source->modelPSF != NULL) continue;
-
         // skip non-astronomical objects (very likely defects)
         if (source->type == PS_SOURCE_DEFECT) continue;
@@ -54 +52 @@
         if (isnan(model->radius)) {
           fprintf (stderr, "error in radius\n");
+          continue;
         }
             

trunk/psphot/src/mark_psf_source.c

@@ -22 +22 @@
     if (source->modelPSF == NULL) return (false);
 
-    if (source->type == PS_SOURCE_PSFSTAR) return (true);
+    // SATSTAR objects are fitted with PSF regardless
     if (source->type == PS_SOURCE_SATSTAR) return (true);
 
@@ -36 +36 @@
     // assign PS_SOURCE_BRIGHTSTAR to bright objects within PSF region of dparams[]
     if ((fabs(nSx) < shapeNsigma) && (fabs(nSy) < shapeNsigma)) {
+      if (source->type == PS_SOURCE_PSFSTAR) return (true);
       source->type = PS_SOURCE_BRIGHTSTAR;
       return (true);
     }
     
+    if (source->type == PS_SOURCE_PSFSTAR) {
+      psLogMsg ("psphot", 3, "PSFSTAR demoted based on dSx, dSy\n");
+    }
+
     if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) {
       source->type = PS_SOURCE_GALAXY;

trunk/psphot/src/psPolynomials.c

@@ -453 +453 @@
                             const psVector* dz)
 {
+    PS_VECTOR_CHECK_NULL(mask, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_NULL(z, NULL);
+    PS_VECTOR_CHECK_NULL(dz, NULL);
+
     psVector *zFit   = NULL;
-    psVector *zResid = NULL;
-    psStats  *stats  = NULL;
+    psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64);
+    psStats  *stats  = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
 
     for (int N = 0; N < 3; N++) {
         poly   = VectorFitPolynomial2DOrd_EAM (poly, mask, x, y, z, dz);
         zFit   = Polynomial2DEvalVectorD (poly, x, y);
-        zResid = (psVector *) psBinaryOp (NULL, (void *) z, "-", (void *) zFit);
-
-        stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+        zResid = (psVector *) psBinaryOp (zResid, (void *) z, "-", (void *) zFit);
+
         stats  = psVectorStats (stats, zResid, NULL, mask, 1);
-        psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g (%d pts)\n", 
-                 stats->clippedMean, stats->clippedStdev, stats->clippedNvalues);
+        psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n", 
+                 stats->sampleMean, stats->sampleStdev);
 
         // set mask if pts are not valid
-        int n = 0;
+        // we are masking out any point which is out of range 
+        // recovery is not allowed with this scheme
         for (int i = 0; i < zResid->n; i++) {
-            if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) {
-              mask->data.U8[i] = 1;
-              continue;
-            } else {
-              mask->data.U8[i] = 0;
-            }         
-            n++;
-        }
-    }
+          if (mask->data.U8[i]) continue;
+          if (fabs(zResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) {
+            mask->data.U8[i] = 1;
+            continue;
+          }       
+        }
+        psFree (zFit);
+    }
+    psFree (zResid);
+    psFree (stats);
     return (poly);
 }

trunk/psphot/src/pspsf.c

@@ -172 +172 @@
     // we are doing a robust fit.  after each pass, we drop points which are 
     // more deviant than three sigma. 
-    // mask is currently updated for each pass. this is inconsistent
+    // mask is currently updated for each pass. this is inconsistent?
 
     for (int i = 0; i < psf->params->n; i++) {

trunk/psphot/src/subtract_psf_source.c

@@ -9 +9 @@
   psSource *source = sources->data[i];
 
+  // non-stellar sources are ignored
   if (source->type == PS_SOURCE_GALAXY) return (false);
   if (source->type == PS_SOURCE_DEFECT) return (false);