Changeset 4115 for trunk/psphot/src

Timestamp:

Jun 4, 2005, 5:22:13 PM (21 years ago)

Author:

eugene

Message:

running sources in brightness order, using a mask

Location:

trunk/psphot/src

Files:

• apply_psf_model.c (modified) (3 diffs)
• by_SN.c (added)
• choose_psf_model.c (modified) (1 diff)
• find_peaks.c (modified) (2 diffs)
• image_stats.c (modified) (1 diff)
• mark_psf_source.c (added)
• psPolynomials.c (modified) (8 diffs)
• psphot.c (modified) (2 diffs)
• psphot.h (modified) (1 diff)
• pspsf.c (modified) (8 diffs)
• subtract_psf_source.c (added)
• subtract_psf_sources.c (modified) (1 diff)
• test_psf_scatter.c (modified) (1 diff)

trunk/psphot/src/apply_psf_model.c

@@ -22 +22 @@
     keep           = psRegionForImage (keep, image, keep);
 
-    // set the object surfact-brightness limit for fitted pixels
+    float shapeNsigma = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA");
+    float snFaint     = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
+
+    // set the object surface-brightness limit for fitted pixels
     float FLUX_LIMIT  = FIT_NSIGMA * sky->sampleStdev;
     psLogMsg ("psphot.apply_psf_model", 3, "fitting pixels with at least %f object counts\n", FLUX_LIMIT);
@@ -31 +34 @@
     for (int i = 0; i < sources->n; i++) {
         psSource *source = sources->data[i];
+
+        // skip the existing PSF stars 
+        // XXX drop this -- refit?
         if (source->modelPSF != NULL) continue;
 
@@ -55 +61 @@
         pmSourceMaskSaturated (source, SATURATE);
 
-        // fit as PSF, not FLT (drop poor fits)
+        // fit as PSF, not FLT (skip poor fits)
         if (!pmSourceFitModel (source, model, true)) continue;
         source->modelPSF = model;
         Niter += model[0].nIter;
         Nfit ++;
+
+        mark_psf_source (source);
+        subtract_psf_source (source);
     }
     psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);

trunk/psphot/src/choose_psf_model.c

@@ -68 +68 @@
     for (int i = 0; i < test->sources->n; i++) {
         psSource *source = test->sources->data[i];
-        source->modelPSF = test->modelPSF->data[i];
+        // drop masked sources from PSFSTAR list
+        if (test->mask->data.U8[i]) {
+          source->type = PS_SOURCE_OTHER;
+          // XXX is this the right type to go to?
+        } else {
+          source->modelPSF = test->modelPSF->data[i];
+        }
     }
     return (test->psf);

trunk/psphot/src/find_peaks.c

@@ -9 +9 @@
 
     // smooth the image 
-    // should we also subtract a super-binned image? (as an option)
+    // should we also subtract a super-binned image? (as an option?)
 
     psTimerStart ("psphot");
@@ -24 +24 @@
 
     psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot"));
-
 
     // find the peaks in the smoothed image 

trunk/psphot/src/image_stats.c

@@ -13 +13 @@
     // pass on the stats
     { 
-        bool   status    = false;
+        bool     status  = false;
+        int      Npixels = image->numRows*image->numCols;
+        int      Nsubset = PS_MIN (Npixels, psMetadataLookupF32 (&status, config, "NSUBSET"));
+        psVector *subset = psVectorAlloc (Nsubset, PS_TYPE_F32);
         psRandom *rnd    = psRandomAlloc (PS_RANDOM_TAUS, 0);
-        int    Nsubset   = psMetadataLookupF32 (&status, config, "NSUBSET");
-        double fSubset   = (double) Nsubset / (image->numRows*image->numCols);
-        psVector *subset = psVectorAlloc (Nsubset, PS_TYPE_F32);
 
-        subset->n = 0;
-        for (int i = 0; i < image->numRows; i++) {
-            for (int j = 0; j < image->numCols; j++) {
-                double frnd = psRandomUniform (rnd);
-                if ((fSubset < 1.0) && (fSubset < frnd)) continue;
-                subset->data.F32[subset->n] = image->data.F32[j][i];
-                subset->n ++;
-                if (subset->n == Nsubset) goto got_subset;
-            }
+        // choose Nsubset points between 0 and Nx*Ny, convert to coords
+        Npixels = image->numRows*image->numCols;
+        for (int i = 0; i < Nsubset; i++) {
+          double frnd = psRandomUniform (rnd);
+          int pixel = Npixels * frnd;
+          int ix = pixel / image->numCols;
+          int iy = pixel % image->numCols;
+          subset->data.F32[i] = image->data.F32[iy][ix];
         }
 
-    got_subset:
         // this should use ROBUST not SAMPLE median
         // robust median is broken in pslib (0.5)

trunk/psphot/src/psPolynomials.c

@@ -191 +191 @@
 // XXX EAM : test version of 1d fitting
 psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly,
-        const psVector* x,
-        const psVector* y,
-        const psVector* yErr)
+                                             psVector *mask,
+                                             const psVector *x,
+                                             const psVector *y,
+                                             const psVector *yErr)
 {
     // I think this is 1 dimension down
@@ -232 +233 @@
     // Build the B and A data structs.
     for (int k = 0; k < x->n; k++) {
+        if ((mask != NULL) && mask->U8[k]) continue;
         xSums = BuildSums1D(xSums, x->data.F64[k], nTerm);
       
@@ -291 +293 @@
 // XXX EAM : test version of 2d fitting
 psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly,
-        const psVector* x,
-        const psVector* y,
-        const psVector* z,
-        const psVector* zErr)
+                                             psVector* mask,
+                                             const psVector* x,
+                                             const psVector* y,
+                                             const psVector* z,
+                                             const psVector* zErr)
 {
     // I think this is 1 dimension down
@@ -322 +325 @@
     // Build the B and A data structs.
     for (int k  = 0; k < x->n; k++) {
+      if ((mask != NULL) && mask->data.U8[i]) continue;
         Sums = BuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm);
       
@@ -391 +395 @@
 }    
 
-psPolynomial2D* RobustFit2D(psPolynomial2D* poly,
+psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly,
                             const psVector* x,
                             const psVector* y,
@@ -413 +417 @@
     for (int N = 0; N < 3; N++) {
         // XXX EAM : this would be better defined with an element mask
-        poly   = VectorFitPolynomial2DOrd_EAM (poly, X, Y, Z, dZ);
+        poly   = VectorFitPolynomial2DOrd_EAM (poly, NULL, X, Y, Z, dZ);
         zFit   = Polynomial2DEvalVectorD (poly, x, y);
         zResid = (psVector *) psBinaryOp (NULL, (void *) z, "-", (void *) zFit);
@@ -424 +428 @@
         int n = 0;
         for (int i = 0; i < zResid->n; i++) {
-            if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) continue;
+            if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) {
+              continue;
+            } 
             X->data.F64[n]  =  x->data.F64[i];
             Y->data.F64[n]  =  y->data.F64[i];
@@ -435 +441 @@
         Z->n = n;
         dZ->n = n;
+    }
+    return (poly);
+}
+
+// XXX EAM : be careful here with F32 vs F64 vectors 
+psPolynomial2D* RobustFit2D(psPolynomial2D* poly,
+                            psVector* mask,
+                            const psVector* x,
+                            const psVector* y,
+                            const psVector* z,
+                            const psVector* dz)
+{
+    psVector *zFit   = NULL;
+    psVector *zResid = NULL;
+    psStats  *stats  = NULL;
+
+    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);
+        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);
+
+        // set mask if pts are not valid
+        int n = 0;
+        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++;
+        }
     }
     return (poly);

trunk/psphot/src/psphot.c

@@ -24 +24 @@
     peaks = find_peaks (image, config, sky);
 
+    // construct sources and measure basic stats
     sources = source_moments (image, config, peaks);
+
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, by_SN);
 
     // use stellar objects SN > PSF_SN_LIM
@@ -35 +39 @@
     // XXX it has some problems
     // test_psf_scatter (sources);
-
-    // identify PSF-like objects (< PSF_SHAPE_NSIGMA from dSx,dSy model)
-    // SN <= FAINT_SN_LIM : FAINTSTAR
-    // SN >  FAINT_SN_LIM : BRIGHTSTAR
-    mark_psf_sources (sources, config);
-
-    // subtract BRIGHTSTAR, PSFSTAR, SATSTAR
-    subtract_psf_sources (sources);
 
     // fit_galaxies (image, config, sources);

trunk/psphot/src/psphot.h

@@ -13 +13 @@
   psArray    *modelFLT;     // model fits, floating parameters 
   psArray    *modelPSF;     // model fits, PSF parameters
+  psVector   *mask;
   psVector   *metric;
   psStats    *metricStats;

trunk/psphot/src/pspsf.c

@@ -36 +36 @@
     test->metricStats = NULL;
     test->metric      = psVectorAlloc (sources->n, PS_TYPE_F64);
+    test->mask        = psVectorAlloc (sources->n, PS_TYPE_U8);
 
     for (int i = 0; i < test->modelFLT->n; i++) {
         test->modelFLT->data[i] = NULL;
         test->modelPSF->data[i] = NULL;
+        test->mask->data.U8[i]  = 0;
     }   
     return (test);
@@ -57 +59 @@
         psModel  *model  = pmSourceModelGuess (source, test->modelType); 
 
-        // fit model as FLT, not PSF (skip poor fits)
-        if (!pmSourceFitModel (source, model, false)) continue;
+        // fit model as FLT, not PSF (mask & skip poor fits)
+        if (!pmSourceFitModel (source, model, false)) {
+          test->mask->data.U8[i] = 1;
+          continue;
+        }
         test->modelFLT->data[i] = model;
         Nflt ++;
@@ -66 +71 @@
 
     // stage 2: construct a psf (pmPSF) from this collection of model fits
-    pmPSFFromModels (test->psf, test->modelFLT);
+    pmPSFFromModels (test->psf, test->modelFLT, test->mask);
 
     // stage 3: refit with fixed shape parameters
@@ -73 +78 @@
         psSource *source = test->sources->data[i];
         psModel  *modelFLT = test->modelFLT->data[i];
-        if (modelFLT == NULL) continue;
+
+        // masked for: bad model fit, outlier in parameters
+        if (test->mask->data.U8[i]) continue;
+        if (modelFLT == NULL) {
+          psLogMsg ("psphot.psftest", 2, "programming error: missing model not masked\n");
+          continue;
+        }
         psModel  *modelPSF = psModelFromPSF (modelFLT, test->psf); // set shape for this model
 
         // fit model as PSF, not FLT (skip poor fits)
-        if (!pmSourceFitModel (source, modelPSF, true)) continue;
+        if (!pmSourceFitModel (source, modelPSF, true)) {
+          test->mask->data.U8[i] = 1;
+          continue;
+        }
         test->modelPSF->data[i] = modelPSF;
         Npsf ++;
@@ -92 +106 @@
         float fitSum = 0;
       
+        // is this metricMask redundant with test->mask ?
         metricMask->data.U8[i] = 1;
         test->metric->data.F64[i] = 0;
 
+        if (test->mask->data.U8[i]) continue;
+
         psModel *model = test->modelPSF->data[i];
-        if (model == NULL) continue;
+        if (model == NULL) {
+          psLogMsg ("psphot.psftest", 2, "programming error: missing model not masked\n");
+          continue;
+        }
 
         psImage *image = ((psSource *)test->sources->data[i])->pixels;
         psImage *mask  = ((psSource *)test->sources->data[i])->mask;
 
+        // this metric is Ap-Fit
         float sky = model->params->data.F32[0];
-
         for (int ix = 0; ix < image->numCols; ix++) {
             for (int iy = 0; iy < image->numRows; iy++) {
@@ -110 +130 @@
             }
         }
-        // this metric is Ap-Fit
         // fprintf (stderr, "%d %f %f  %f\n", i, obsSum, fitSum, test->metric->data.F64[i]);
+        // keep the good metrics
         if ((fitSum > 0) && (obsSum > 0)) {
             metricMask->data.U8[i] = 0;
             test->metric->data.F64[i] = -2.5*log10(obsSum/fitSum);
-        }       
+        } else {
+          test->mask->data.U8[i] = 1;
+        }
     }
 
@@ -128 +150 @@
 // input: an array of psModels, pre-allocated psf
 // some of the array entries may be NULL, ignore them
-bool pmPSFFromModels (pmPSF *psf, psArray *models) {
+bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask) {
 
     int n;
 
-    // construct the x and y vectors from the collection of objects
-    // count the number of valid model measurements
+    // construct the fit vectors from the collection of objects
+    // use the mask to ignore missing fits 
     psVector *x = psVectorAlloc (models->n, PS_TYPE_F64);
     psVector *y = psVectorAlloc (models->n, PS_TYPE_F64);
-    n = 0;
+    psVector *z = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64);
+
     for (int i = 0; i < models->n; i++) {
         psModel *model = models->data[i];
@@ -142 +166 @@
 
         // XXX EAM : this is fragile: x and y must be stored consistently at 2,3
-        x->data.F64[n] = model->params->data.F32[2];
-        y->data.F64[n] = model->params->data.F32[3];
-        n++;
-    }
-    x->n = y->n = n;
+        x->data.F64[i] = model->params->data.F32[2];
+        y->data.F64[i] = model->params->data.F32[3];
+    }
 
     // we are doing a robust fit.  after each pass, we drop points which are 
-    // more deviant than three sigma.
-    psVector *z = psVectorAlloc (x->n, PS_TYPE_F64);
-    psVector *dz = psVectorAlloc (x->n, PS_TYPE_F64);
+    // more deviant than three sigma. 
+    // mask is currently updated for each pass. this is inconsistent
+
     for (int i = 0; i < psf->params->n; i++) {
-        n = 0;
         for (int j = 0; j < models->n; j++) {
             psModel *model = models->data[j];
             if (model == NULL) continue;
-            z->data.F64[n] = model->params->data.F32[i + 4];
-            dz->data.F64[n] = 1;
-            n++;
+            z->data.F64[j] = model->params->data.F32[i + 4];
+            dz->data.F64[j] = 1;
             // XXX EAM : need to use actual errors?
             // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4
         }
-        if (n != x->n) psAbort ("pmPSFFromModels", "programming error: mismatched number of NULL models\n");
-        z->n = n;
-        dz->n = n;
-
-        psf->params->data[i] = RobustFit2D (psf->params->data[i], x, y, z, dz);
+        // if (n != x->n) psAbort ("pmPSFFromModels", "programming error: mismatched number of NULL models\n");
+
+        psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
         // psPolynomial2DDump (psf->params->data[i]);
     }

trunk/psphot/src/subtract_psf_sources.c

@@ -6 +6 @@
 bool subtract_psf_sources (psArray *sources)
 { 
+
+    sources = psArraySort (sources, by_SN);
+
     for (int i = 0; i < sources->n; i++) {
         float sky;

trunk/psphot/src/test_psf_scatter.c

@@ -71 +71 @@
 
     psPolynomial1D *dsxLine = psPolynomial1DAlloc (2, PS_POLYNOMIAL_ORD);
-    psVectorFitPolynomial1D (dsxLine, snBin, dsxBin, NULL);
+    psVectorFitPolynomial1D (dsxLine, NULL, snBin, dsxBin, NULL);
     psPolynomial1D *dsyLine = psPolynomial1DAlloc (2, PS_POLYNOMIAL_ORD);
-    psVectorFitPolynomial1D (dsyLine, snBin, dsyBin, NULL);
+    psVectorFitPolynomial1D (dsyLine, NULL, snBin, dsyBin, NULL);
     // these should have a slope of 1.0
     psPolynomial1DDump (dsxLine);