Changeset 33415 for branches/meh_branches/ppstack_test/psphot/src/psphotExtendedSourceAnalysis.c

Timestamp:

Mar 5, 2012, 5:19:48 PM (14 years ago)

Author:

mhuber

Message:

merging latest r33407 trunk changes to meh_branches/ppstack_test

Location:

branches/meh_branches/ppstack_test

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (2 props)
• psphot/src/psphotExtendedSourceAnalysis.c (modified) (8 diffs)

branches/meh_branches/ppstack_test/psphot/src

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest

branches/meh_branches/ppstack_test/psphot/src/psphotExtendedSourceAnalysis.c

@@ -2 +2 @@
 
 // measure the elliptical radial profile and use this to measure the petrosian parameters for the sources
-// XXX this function needs to be threaded
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -8 +7 @@
 {
     bool status = true;
+
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Extended Source Analysis (Petrosians) ---");
 
     // select the appropriate recipe information
@@ -59 +61 @@
     }
 
-    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
-    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
-    assert (maskVal);
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
 
     // get the sky noise from the background analysis; if this is missing, get the user-supplied value
@@ -70 +74 @@
     }
 
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, pmSourceSortByFlux);
+
+    // option to limit analysis to a specific region
+    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
+    psRegion *AnalysisRegion = psRegionAlloc(0,0,0,0);
+    *AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
+    if (psRegionIsNaN (*AnalysisRegion)) psAbort("analysis region mis-defined");
+
+    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
+    int Cx = 1, Cy = 1;
+    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
+
+    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
+
+    for (int i = 0; i < cellGroups->n; i++) {
+
+        psArray *cells = cellGroups->data[i];
+
+        for (int j = 0; j < cells->n; j++) {
+
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_EXTENDED_ANALYSIS");
+
+            psArrayAdd(job->args, 1, readout);
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, AnalysisRegion);
+            psArrayAdd(job->args, 1, recipe);
+
+            PS_ARRAY_ADD_SCALAR(job->args, skynoise, PS_TYPE_F32);
+
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npetro
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nannuli
+
+// 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 (!psphotExtendedSourceAnalysis_Threaded(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree(AnalysisRegion);
+                return false;
+            }
+            psScalar *scalar = NULL;
+            scalar = job->args->data[5];
+            Next += scalar->data.S32;
+            scalar = job->args->data[6];
+            Npetro += scalar->data.S32;
+            scalar = job->args->data[7];
+            Nannuli += scalar->data.S32;
+            psFree(job);
+# endif
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false, true)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            psFree(AnalysisRegion);
+            return false;
+        }
+
+        // we have only supplied one type of job, so we can assume the types here
+        psThreadJob *job = NULL;
+        while ((job = psThreadJobGetDone()) != NULL) {
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            } else {
+                psScalar *scalar = NULL;
+                scalar = job->args->data[5];
+                Next += scalar->data.S32;
+                scalar = job->args->data[6];
+                Npetro += scalar->data.S32;
+                scalar = job->args->data[7];
+                Nannuli += scalar->data.S32;
+            }
+            psFree(job);
+        }
+    }
+    psFree (cellGroups);
+    psFree(AnalysisRegion);
+
+    psLogMsg ("psphot", PS_LOG_WARN, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
+    psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
+    psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
+
+    psphotVisualShowResidualImage (readout, false);
+
+    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
+    if (doPetrosian) {
+        psphotVisualShowPetrosians (sources);
+    }
+
+    return true;
+}
+
+bool psphotExtendedSourceAnalysis_Threaded (psThreadJob *job) {
+
+    bool status;
+
+    int Next = 0;
+    int Npetro = 0;
+    int Nannuli = 0;
+
+    // arguments: readout, sources, models, region, psfSize, maskVal, markVal
+    pmReadout *readout      = job->args->data[0];
+    psArray *sources        = job->args->data[1];
+    psRegion *region        = job->args->data[2];
+    psMetadata *recipe      = job->args->data[3];
+
+    float skynoise          = PS_SCALAR_VALUE(job->args->data[4],F32);
+
     // S/N limit to perform full non-linear fits
     float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
@@ -78 +198 @@
     bool doPetroStars   = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS");
 
-    // source analysis is done in S/N order (brightest first)
-    sources = psArraySort (sources, pmSourceSortByFlux);
-
-    // option to limit analysis to a specific region
-    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
-    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
-    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    assert (maskVal);
 
     // choose the sources of interest
@@ -90 +206 @@
 
         pmSource *source = sources->data[i];
+
+        // if we have checked the source validity on the basis of the object set, then 
+        // we either skip these tests below or we skip the source completely
+        if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_SKIP) continue;
+        if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_KEEP) goto keepSource;
 
         // skip PSF-like and non-astronomical objects
@@ -104 +225 @@
 
         // limit selection to some SN limit
-        assert (source->peak); // how can a source not have a peak?
-        if (sqrt(source->peak->detValue) < SN_LIM) continue;
-
-        // limit selection by analysis region
-        if (source->peak->x < AnalysisRegion.x0) continue;
-        if (source->peak->y < AnalysisRegion.y0) continue;
-        if (source->peak->x > AnalysisRegion.x1) continue;
-        if (source->peak->y > AnalysisRegion.y1) continue;
+        // assert (source->peak); // how can a source not have a peak?
+        // limit selection to some SN limit
+        bool skipSource = false;
+        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+            skipSource = (source->moments->KronFlux < SN_LIM * source->moments->KronFluxErr);
+        } else {
+            skipSource = (sqrt(source->peak->detValue) < SN_LIM);
+        }
+        if (skipSource) continue;
+
+        // limit selection by analysis region (this automatically apply
+        if (source->peak->x < region->x0) continue;
+        if (source->peak->y < region->y0) continue;
+        if (source->peak->x > region->x1) continue;
+        if (source->peak->y > region->y1) continue;
+
+    keepSource:
 
         // replace object in image
@@ -159 +289 @@
     }
 
-    psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
-    psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
-    psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
-
-    psphotVisualShowResidualImage (readout, false);
-
-    if (doPetrosian) {
-        psphotVisualShowPetrosians (sources);
-    }
-
-    // fprintf (stderr, "xsrc : tried %ld objects\n", sources->n);
+    psScalar *scalar = NULL;
+
+    // change the value of a scalar on the array (wrap this and put it in psArray.h)
+    scalar = job->args->data[5];
+    scalar->data.S32 = Next;
+
+    scalar = job->args->data[6];
+    scalar->data.S32 = Npetro;
+
+    scalar = job->args->data[7];
+    scalar->data.S32 = Nannuli;
 
     return true;

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest