Changeset 32023 for branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceAnalysis.c

Timestamp:

Aug 3, 2011, 6:29:49 PM (15 years ago)

Author:

eugene

Message:

thread psphotExtendedAnalsis and psphotRadialApertures; subtract the fake source image so output residuals image is clean

File:

• branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceAnalysis.c (modified) (6 diffs)

branches/eam_branches/ipp-20110710/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
@@ -59 +58 @@
     }
 
-    // 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 +71 @@
     }
 
+    // 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)) {
+            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_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);
+
+    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 +195 @@
     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
@@ -108 +221 @@
 
         // 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;
+        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;
 
         // replace object in image
@@ -159 +272 @@
     }
 
-    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;