Changeset 32023 for branches/eam_branches/ipp-20110710/psphot/src/psphotRadialApertures.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/psphotRadialApertures.c (modified) (5 diffs)

branches/eam_branches/ipp-20110710/psphot/src/psphotRadialApertures.c

@@ -63 +63 @@
     }
 
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
+
+    // source analysis is done in S/N order (brightest first)
+    // XXX are we getting the objects out of order? does it matter?
+    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_RADIAL_APERTURES");
+
+            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, 0, PS_TYPE_S32); // this is used as a return value for Nradial
+
+// 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 (!psphotRadialApertures_Threaded(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree(AnalysisRegion);
+                return false;
+            }
+            psScalar *scalar = NULL;
+            scalar = job->args->data[4];
+            Nradial += 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[4];
+                Nradial += scalar->data.S32;
+            }
+            psFree(job);
+        }
+    }
+    psFree (cellGroups);
+    psFree(AnalysisRegion);
+
+    psLogMsg ("psphot", PS_LOG_INFO, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);
+    return true;
+}
+ 
+bool psphotRadialApertures_Threaded (psThreadJob *job) {
+
+    bool status;
+    int Nradial = 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];
+
     // radMax stores the upper bounds of the annuli
     // XXX keep the same name here as for the petrosian / elliptical apertures?
@@ -68 +164 @@
     psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
     psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
-    float outerRadius = radMax->data.F32[radMax->n - 1];
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
@@ -77 +172 @@
     float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
 
-    // source analysis is done in S/N order (brightest first)
-    // XXX are we getting the objects out of order? does it matter?
-    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");
+    float outerRadius = radMax->data.F32[radMax->n - 1];
 
     // choose the sources of interest
@@ -104 +192 @@
 
         // 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;
 
         // allocate pmSourceExtendedParameters, if not already defined
@@ -145 +233 @@
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     }
-
-    psLogMsg ("psphot", PS_LOG_INFO, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);
+    psScalar *scalar = job->args->data[4];
+    scalar->data.S32 = Nradial;
+
     return true;
 }