Changeset 35769 for trunk/psphot/src/psphotModelBackground.c

Timestamp:

Jul 3, 2013, 2:43:13 PM (13 years ago)

Author:

eugene

Message:

thread psphotModelBackground, psphotAddOrSubNoise, psphotGuessModel; add recipe options to control selection of extended models fitted; iterate on the trail window a bit

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/psphotModelBackground.c (modified) (9 diffs)

trunk/psphot/src/psphotModelBackground.c

@@ -28 +28 @@
 }
 
+// track background cell failures for log reporting purposes
+static int nFailures = 0;
 
 // Generate the background model
@@ -33 +35 @@
 // linear interpolation to generate full-scale model
 //
-// NOTE that the 'analysis' metedata pass in here is used to store the binning information.
+// NOTE that the 'analysis' metedata passed in here is used to store the binning information.
 // This may be the analysis for this readout, but it may be the analysis for the pmFPAfile
 // corresponding to the model.  Other information about the background model is saved on the
@@ -109 +111 @@
 
     const psStatsOptions statsOption = statsOptionLocation | statsOptionWidth;
-    psStats *stats = psStatsAlloc (statsOption);
+    psStats *statsDefaults = psStatsAlloc (statsOption);
 
     // set range for old-version of sky statistic
     if (statsOptionLocation & PS_STAT_ROBUST_QUARTILE) {
-        stats->min = 0.25;
-        stats->max = 0.75;
+        statsDefaults->min = 0.25;
+        statsDefaults->max = 0.75;
     }
 
     // set user-option for number of pixels per region
-    stats->nSubsample = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
+    statsDefaults->nSubsample = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
     if (!status) {
-        stats->nSubsample = 1000;
+        statsDefaults->nSubsample = 1000;
     }
 
     // optionally set the binsize
-    stats->binsize = psMetadataLookupF32 (&status, recipe, "SKY_HISTOGRAM_BINS");
+    statsDefaults->binsize = psMetadataLookupF32 (&status, recipe, "SKY_HISTOGRAM_BINS");
     if (status) {
-        stats->options |= PS_STAT_USE_BINSIZE;
+        statsDefaults->options |= PS_STAT_USE_BINSIZE;
     }
 
     // optionally set the sigma clipping
-    stats->clipSigma = psMetadataLookupF32 (&status, recipe, "SKY_CLIP_SIGMA");
+    statsDefaults->clipSigma = psMetadataLookupF32 (&status, recipe, "SKY_CLIP_SIGMA");
     if (!status) {
-        if (stats->options & PS_STAT_FITTED_MEAN) {
-            stats->clipSigma = 1.0;
+        if (statsDefaults->options & PS_STAT_FITTED_MEAN) {
+            statsDefaults->clipSigma = 1.0;
         } else {
-            stats->clipSigma = 3.0;
+            statsDefaults->clipSigma = 3.0;
         }
     }
-
-    // stats is not initialized by psStats???  use this to save the input options
-    // XXX re-work this to use the new psStatsInit function
-    psStats *statsDefaults = psStatsAlloc (statsOption);
-    *statsDefaults = *stats;
 
     // we save the binning structure for use in psphotMagnitudes
@@ -154 +151 @@
 
     // XXXX we can thread this here by running blocks in parallel
-
-    int nFailures = 0;
     psImageBackgroundInit();
+    nFailures = 0; // reset for this pass
 
     // we have Nx * Ny model points, but we can use a window which is larger (or smaller) than
@@ -163 +159 @@
 
     // measure clipped median for subimages
-    psRegion ruffRegion = psRegionSet (0,0,0,0);
-    psRegion fineRegion = psRegionSet (0,0,0,0);
     for (int iy = 0; iy < model->numRows; iy++) {
         for (int ix = 0; ix < model->numCols; ix++) {
-
-            // convert the ruff grid cell to the equivalent fine grid cell
-            ruffRegion = psRegionSet (ix + 0.5 - 0.5*dXsample, ix + 0.5 + 0.5*dXsample, 
-                                      iy + 0.5 - 0.5*dYsample, iy + 0.5 + 0.5*dYsample);
-            fineRegion = psImageBinningSetFineRegion (binning, ruffRegion);
-            fineRegion = psRegionForImage (image, fineRegion);
-
-            psImage *subset  = psImageSubset (image, fineRegion);
-            if (!subset->numCols || !subset->numRows) {
-                psFree (subset);
-                continue;
+            
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_MODEL_BACKGROUND");
+
+            psArrayAdd(job->args, 1, image);
+            psArrayAdd(job->args, 1, mask);
+            psArrayAdd(job->args, 1, binning);
+            psArrayAdd(job->args, 1, rng);
+            psArrayAdd(job->args, 1, statsDefaults);
+
+            psArrayAdd(job->args, 1, modelData);
+            psArrayAdd(job->args, 1, modelStdevData);
+
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+
+            PS_ARRAY_ADD_SCALAR(job->args, ix, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, iy, PS_TYPE_S32);
+
+            PS_ARRAY_ADD_SCALAR(job->args, dXsample, PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, dYsample, PS_TYPE_F32);
+
+            PS_ARRAY_ADD_SCALAR(job->args, statsOptionLocation, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, statsOptionWidth, PS_TYPE_S32);
+
+            PS_ARRAY_ADD_SCALAR(job->args, NAN, PS_TYPE_F32); // this is used as a return value for dQvalue
+
+# if (1)
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                return NULL;
             }
-            psImage *submask = psImageSubset (mask, fineRegion);
-
-            // reset the default values
-            *stats = *statsDefaults;
-
-            // Use the selected background statistic for the first pass
-            // If it fails, fall back on the "ROBUST_MEDIAN" version
-            // If both fail, set the pixel to NAN and (below) interpolate
-            // XXX psImageBackground will probably be renamed psImageStats
-            // XXX don't bother trying if there are no valid pixels...
-
-            psVector *sample = NULL;
-
-            // turn on stats tracing in desired cells
-            # if (0)
-            psMetadata *plots = psMetadataLookupPtr (&status, recipe, "DIAGNOSTIC.PLOTS");
-            assert (plots);
-
-            int xPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.X");
-            assert (status);
-            int yPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.Y");
-            assert (status);
-
-            bool gotX = (xPlot < 0) || (xPlot == ix);
-            bool gotY = (yPlot < 0) || (yPlot == iy);
-
-            if (gotX && gotY) {
-                (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 6);
-                (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 6);
-            } else {
-                (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 0);
-                (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 0);
+# else
+            if (!psphotModelBackground_Threaded(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Failure to model background.");
+                return NULL;
             }
-            # endif
-
-            if (psImageBackground(stats, &sample, subset, submask, maskVal, rng)) {
-                if (stats->options & PS_STAT_ROBUST_QUARTILE) {
-                    modelData[iy][ix] = stats->robustMedian;
-                } else {
-                    modelData[iy][ix] = psStatsGetValue(stats, statsOptionLocation);
-                }
-                modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth);
-                // fprintf (stderr, "dQ : %f - %f - %f = %f\n", stats->robustLQ, stats->robustMedian, stats->robustUQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian);
-                psVectorAppend (dQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian);
-
-                // supply sample to plotting routing
-                psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample);
-            } else {
-                psStatsOptions currentOptions = stats->options;
-                stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV;
-                if (!psImageBackground(stats, &sample, subset, submask, maskVal, rng)) {
-                    if ((nFailures < 3) || (nFailures % 100 == 0)) {
-                        psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for "
-                                  "(%dx%d, (row0,col0) = (%d,%d)",
-                                  subset->numRows, subset->numCols, subset->row0, subset->col0);
-                    }
-                    nFailures ++;
-                    modelData[iy][ix] = modelStdevData[iy][ix] = NAN;
-                } else {
-                    modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN);
-                    modelStdevData[iy][ix] = psStatsGetValue(stats, PS_STAT_ROBUST_STDEV);
-
-                    // supply sample to plotting routing
-                    psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample);
-                }
-                // drop errors caused by psImageBackground failures
-                // XXX we probably should trap and exit on serious failures
-                psErrorClear();
-                stats->options = currentOptions;
-            }
-            psFree(sample);
-            modelData[iy][ix] += SKY_BIAS;
-            psFree (subset);
-            psFree (submask);
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            } else {
+                psScalar *scalar = job->args->data[14];
+                float dQvalue = scalar->data.F32;
+                psVectorAppend (dQ, dQvalue);
+            }
+            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.");
+        return NULL;
+    }
+
+    // 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 = job->args->data[14];
+            float dQvalue = scalar->data.F32;
+            psVectorAppend (dQ, dQvalue);
+        }
+        psFree(job);
     }
 
@@ -302 +274 @@
     if (Count == 0) {
         psError (PSPHOT_ERR_DATA, true, "failed to build background image");
-        psFree(stats);
         psFree(statsDefaults);
         psFree(binning);
@@ -318 +289 @@
             modelData[iy][ix] = Value;
             modelStdevData[iy][ix] = ValueStdev;
+        }
+    }
+
+    // apply artificial offset if desired
+    for (int iy = 0; iy < model->numRows; iy++) {
+        for (int ix = 0; ix < model->numCols; ix++) {
+            modelData[iy][ix] += SKY_BIAS;
         }
     }
@@ -364 +342 @@
     psFree(dQ);
 
-    psFree(stats);
     psFree(statsDefaults);
     psFree(binning);
@@ -433 +410 @@
     return true;
 }
+
+bool psphotModelBackground_Threaded (psThreadJob *job) {
+
+    psImage *image          = job->args->data[0];
+    psImage *mask           = job->args->data[1];
+
+    psImageBinning *binning = job->args->data[2];
+
+    psRandom *rng           = job->args->data[3];
+    psStats *statsDefaults  = job->args->data[4];
+
+    psF32 **modelData       = job->args->data[5];
+    psF32 **modelStdevData  = job->args->data[6];
+
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA);
+
+    int ix                  = PS_SCALAR_VALUE(job->args->data[8],S32);
+    int iy                  = PS_SCALAR_VALUE(job->args->data[9],S32);
+
+    float dXsample          = PS_SCALAR_VALUE(job->args->data[10],F32);
+    float dYsample          = PS_SCALAR_VALUE(job->args->data[11],F32);
+
+    psStatsOptions statsOptionLocation = PS_SCALAR_VALUE(job->args->data[12],S32);
+    psStatsOptions statsOptionWidth    = PS_SCALAR_VALUE(job->args->data[13],S32);
+
+    // convert the ruff grid cell to the equivalent fine grid cell
+    psRegion ruffRegion = psRegionSet (ix + 0.5 - 0.5*dXsample, ix + 0.5 + 0.5*dXsample, iy + 0.5 - 0.5*dYsample, iy + 0.5 + 0.5*dYsample);
+    psRegion fineRegion = psImageBinningSetFineRegion (binning, ruffRegion);
+    fineRegion = psRegionForImage (image, fineRegion);
+
+    psImage *subset  = psImageSubset (image, fineRegion);
+    if (!subset->numCols || !subset->numRows) {
+        psFree (subset);
+        return false; // XXX do we / should we fail on this?
+    }
+    psImage *submask = psImageSubset (mask, fineRegion);
+
+    psStats *stats = psStatsAlloc (PS_STAT_NONE);
+    *stats = *statsDefaults;
+
+    // Use the selected background statistic for the first pass
+    // If it fails, fall back on the "ROBUST_MEDIAN" version
+    // If both fail, set the pixel to NAN and (later) interpolate
+
+    psVector *sample = NULL;
+    float dQvalue = NAN;
+
+    if (psImageBackground(stats, &sample, subset, submask, maskVal, rng)) {
+        if (stats->options & PS_STAT_ROBUST_QUARTILE) {
+            modelData[iy][ix] = stats->robustMedian;
+        } else {
+            modelData[iy][ix] = psStatsGetValue(stats, statsOptionLocation);
+        }
+        modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth);
+
+        // fprintf (stderr, "dQ : %f - %f - %f = %f\n", stats->robustLQ, stats->robustMedian, stats->robustUQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian);
+        // XXX this operation is not thread safe -- move out somehow...
+        // psVectorAppend (dQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian);
+        dQvalue = stats->robustUQ + stats->robustLQ - 2*stats->robustMedian;
+        // return dQvalue to main thread
+
+        // supply sample to plotting routing
+        // only allow in a non-threaded context -- can plot more than one cell
+        // psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample);
+    } else {
+        // psStatsOptions currentOptions = stats->options;
+        stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV;
+        if (!psImageBackground(stats, &sample, subset, submask, maskVal, rng)) {
+            if ((nFailures < 3) || (nFailures % 100 == 0)) {
+                psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for "
+                          "(%dx%d, (row0,col0) = (%d,%d)",
+                          subset->numRows, subset->numCols, subset->row0, subset->col0);
+            }
+            nFailures ++; // static
+            modelData[iy][ix] = modelStdevData[iy][ix] = NAN;
+        } else {
+            modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN);
+            modelStdevData[iy][ix] = psStatsGetValue(stats, PS_STAT_ROBUST_STDEV);
+
+            // supply sample to plotting routing
+            // only allow in a non-threaded context -- can plot more than one cell
+            // psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample);
+        }
+        // drop errors caused by psImageBackground failures
+        // NOTE : psStats raises errors when it cannot find a solution; this 
+        // probably should not be errors but exit stats info only, but c'est la code
+        // XXX we probably should trap and exit on serious failures
+        psErrorClear();
+        // stats->options = currentOptions; // this is not needed (set by init above)
+    }
+    psFree (stats);
+    psFree (sample);
+    psFree (subset);
+    psFree (submask);
+
+    // return the dQvalue to the calling thread
+    psScalar *scalar = job->args->data[14];
+    scalar->data.F32 = dQvalue;
+
+    return true;
+}
+
+// code for in-line plotting from above
+// turn on stats tracing in desired cells
+// XXX this code may need to be re-worked for a threaded context
+# if (0)
+    psMetadata *plots = psMetadataLookupPtr (&status, recipe, "DIAGNOSTIC.PLOTS");
+    assert (plots);
+
+    int xPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.X");
+    assert (status);
+    int yPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.Y");
+    assert (status);
+
+    bool gotX = (xPlot < 0) || (xPlot == ix);
+    bool gotY = (yPlot < 0) || (yPlot == iy);
+
+    if (gotX && gotY) {
+        (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 6);
+        (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 6);
+    } else {
+        (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 0);
+        (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 0);
+    }
+# endif