Changeset 31154 for trunk/psphot/src/psphotExtendedSourceFits.c

Timestamp:

Apr 4, 2011, 1:12:39 PM (15 years ago)

Author:

eugene

Message:

use the smoothed image for footprint culling; use a more stringent culling for saturated stars; more tweaks to get sat stars to be fitted; various updates to psphotStack; unify psphotImageLoop varients; be a bit careful about number of stars in a PSF clump

Location:

trunk/psphot

Files:

• . (modified) (2 props)
• src/psphotExtendedSourceFits.c (modified) (12 diffs)

trunk/psphot

@@ -19 +19 @@
 psphot-config
 Doxyfile
+a.out.dSYM

@@ -19 +19 @@
 psphot-config
 Doxyfile
+a.out.dSYM

trunk/psphot/src/psphotExtendedSourceFits.c

@@ -43 +43 @@
     int NplainPass = 0;
     int Nfaint = 0;
+    int Nfail = 0;
 
     psTimerStart ("psphot.extended");
@@ -145 +146 @@
 
     // source analysis is done in S/N order (brightest first)
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
@@ -176 +177 @@
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nplain
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NplainPass
-            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfain
-
-            if (false && !psThreadJobAddPending(job)) {
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfaint
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
+
+// set this to 0 to run without threading
+# if (1)            
+            if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
                 psFree(AnalysisRegion);
                 return false;
-            } else {
-                if (!psphotExtendedSourceFits_Threaded(job)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                    psFree(AnalysisRegion);
-                    return false;
-                }
-                psScalar *scalar = NULL;
-                scalar = job->args->data[7];
-                Next += scalar->data.S32;
-                scalar = job->args->data[8];
-                Nconvolve += scalar->data.S32;
-                scalar = job->args->data[9];
-                NconvolvePass += scalar->data.S32;
-                scalar = job->args->data[10];
-                Nplain += scalar->data.S32;
-                scalar = job->args->data[11];
-                NplainPass += scalar->data.S32;
-                scalar = job->args->data[12];
-                Nfaint += scalar->data.S32;
-                psFree(job);
+            } 
+# else
+            if (!psphotExtendedSourceFits_Threaded(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree(AnalysisRegion);
+                return false;
             }
-        }
+            psScalar *scalar = NULL;
+            scalar = job->args->data[7];
+            Next += scalar->data.S32;
+            scalar = job->args->data[8];
+            Nconvolve += scalar->data.S32;
+            scalar = job->args->data[9];
+            NconvolvePass += scalar->data.S32;
+            scalar = job->args->data[10];
+            Nplain += scalar->data.S32;
+            scalar = job->args->data[11];
+            NplainPass += scalar->data.S32;
+            scalar = job->args->data[12];
+            Nfaint += scalar->data.S32;
+            scalar = job->args->data[13];
+            Nfail += scalar->data.S32;
+            psFree(job);
+# endif
+        }
 
         // wait for the threads to finish and manage results
@@ -231 +238 @@
                 scalar = job->args->data[12];
                 Nfaint += scalar->data.S32;
+                scalar = job->args->data[13];
+                Nfail += scalar->data.S32;
             }
             psFree(job);
@@ -238 +247 @@
     psFree(AnalysisRegion);
 
-    psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
+    psLogMsg ("psphot", PS_LOG_INFO, "extended source model fits: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d convolved models (%d passed)\n", Nconvolve, NconvolvePass);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models (%d passed)\n", Nplain, NplainPass);
-    psLogMsg ("psphot", PS_LOG_INFO, "  %d too faint to fit\n", Nfaint);
+    psLogMsg ("psphot", PS_LOG_INFO, "  %d too faint to fit, %d failed\n", Nfaint, Nfail);
     return true;
 }
@@ -250 +259 @@
     bool status;
     int Next = 0;
+    int Nfaint = 0;
+    int Nfail = 0;
     int Nconvolve = 0;
     int NconvolvePass = 0;
     int Nplain = 0;
-    int Nfaint = 0;
     int NplainPass = 0;
     bool savePics = false;
@@ -268 +278 @@
     psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
 
+    // pthread_t tid = pthread_self();     // Thread identifier
+
     // Define source fitting parameters for extended source fits
     pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
-    fitOptions->mode          = PM_SOURCE_FIT_EXT;
+    fitOptions->mode           = PM_SOURCE_FIT_EXT;
+    fitOptions->saveCovariance = true;  // XXX make this a user option?
+    fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
+
     // XXX for now, use the defaults for the rest:
     // fitOptions->nIter         = fitIter;
@@ -296 +311 @@
         // XXX use the parameters guessed from moments
         // if (source->modelEXT == NULL) continue;
+
+        // fprintf (stderr, "fit %d,%d in thread %d\n", source->peak->x, source->peak->y, (int) tid);
 
         // replace object in image
@@ -366 +383 @@
           // limit selection to some SN limit
           assert (source->peak); // how can a source not have a peak?
-          if (source->peak->SN < SNlim) {
+          if (sqrt(source->peak->detValue) < SNlim) {
               Nfaint ++;
               continue;
@@ -388 +405 @@
               if (!modelFit) {
                   psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->Mx, source->moments->My);
+                  Nfail ++;
                   continue;
               }
@@ -402 +420 @@
               if (!modelFit) {
                   psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My);
+                  Nfail ++;
                   continue;
               }
@@ -519 +538 @@
     scalar->data.S32 = Nfaint;
 
+    scalar = job->args->data[13];
+    scalar->data.S32 = Nfail;
+
     return true;
 }

@@ -19 +19 @@
 psphot-config
 Doxyfile
+a.out.dSYM