Changeset 21166

Timestamp:

Jan 24, 2009, 10:54:29 AM (17 years ago)

Author:

eugene

Message:

clean up threading model for psphotGuessModel (defined standard way to split sources into regions); extending the threading to psphotMagnitudes, psphotApReset, psphotBlendFit, psphotSourceStats

Location:

trunk/psphot/src

Files:

• Makefile.am (modified) (1 diff)
• psphot.h (modified) (3 diffs)
• psphotApResid.c (modified) (7 diffs)
• psphotArguments.c (modified) (1 diff)
• psphotBlendFit.c (modified) (8 diffs)
• psphotCleanup.c (modified) (1 diff)
• psphotGuessModels.c (modified) (7 diffs)
• psphotMagnitudes.c (modified) (4 diffs)
• psphotReadout.c (modified) (5 diffs)
• psphotReadoutFindPSF.c (modified) (1 diff)
• psphotSetThreads.c (modified) (2 diffs)
• psphotSourceFits.c (modified) (3 diffs)
• psphotSourceStats.c (modified) (7 diffs)
• psphotThreadTools.c (added)

trunk/psphot/src/Makefile.am

@@ -79 +79 @@
         psphotMakeFluxScale.c          \
         psphotCheckStarDistribution.c  \
+        psphotThreadTools.c            \
         psphotAddNoise.c
 

trunk/psphot/src/psphot.h

@@ -12 +12 @@
 
 #define PSPHOT_RECIPE_PSF_FAKE_ALLOW "PSF.FAKE.ALLOW" // Name for recipe component permitting fake PSFs
-
-typedef struct {
-    pmReadout *readout;
-    psArray *sources;
-    pmPSF *psf;
-    psRegion *region;
-    psMaskType maskVal;
-    psMaskType markVal;
-} psphotGuessModelForRegionArgs;
 
 // top-level psphot functions
@@ -43 +34 @@
 pmDetections   *psphotFindDetections (pmDetections *detections, pmReadout *readout, psMetadata *recipe);
 
-psArray        *psphotSourceStats (pmReadout *readout, psMetadata *recipe, pmDetections *detections);
 bool            psphotRoughClass (pmReadout *readout, psArray *sources, psMetadata *recipe, const bool findPsfClump);
 bool            psphotBasicDeblend (psArray *sources, psMetadata *recipe);
@@ -50 +40 @@
 bool            psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources);
 bool            psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final);
-bool            psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf);
-bool            psphotGuessModelForRegion (psphotGuessModelForRegionArgs *args);
-bool            psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf);
 bool            psphotReplaceUnfitSources (psArray *sources, psMaskType maskVal);
 bool            psphotReplaceAllSources (psArray *sources, psMetadata *recipe);
 bool            psphotRemoveAllSources (psArray *sources, psMetadata *recipe);
-bool            psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf);
-bool            psphotMagnitudes (pmConfig *config, const pmFPAview *view, psArray *sources, psMetadata *recipe, pmPSF *psf);
+
+bool            psphotBlendFit (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf);
+bool            psphotBlendFit_Threaded (psThreadJob *job);
+
+psArray        *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections);
+bool            psphotSourceStats_Threaded (psThreadJob *job);
+
+bool            psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf);
+bool            psphotGuessModel_Threaded (psThreadJob *job);
+
+bool            psphotMagnitudes (pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf);
+bool            psphotMagnitudes_Threaded (psThreadJob *job);
+
+bool            psphotApResid (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf);
+bool            psphotApResidMags_Threaded (psThreadJob *job);
+
 bool            psphotSkyReplace (pmConfig *config, const pmFPAview *view);
 bool            psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe);
 bool            psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe);
+
+// thread-related:
 bool            psphotSetThreads ();
+bool            psphotChooseCellSizes (int *Cx, int *Cy, pmReadout *readout, int nThreads);
+bool            psphotCoordToCell (int *group, int *cell, float x, float y, int Cx, int Cy);
+psArray        *psphotAssignSources (int Cx, int Cy, psArray *sources);
 
 // used by psphotFindDetections

trunk/psphot/src/psphotApResid.c

@@ -4 +4 @@
 // measure the aperture residual statistics and 2D variations
 
-bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf)
+bool psphotApResid (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf)
 {
     int Nfail = 0;
@@ -13 +13 @@
     pmSource *source;
 
-    PS_ASSERT_PTR_NON_NULL(psf, false);
+    PS_ASSERT_PTR_NON_NULL(config, false);
     PS_ASSERT_PTR_NON_NULL(readout, false);
     PS_ASSERT_PTR_NON_NULL(sources, false);
-    PS_ASSERT_PTR_NON_NULL(recipe, false);
+    PS_ASSERT_PTR_NON_NULL(psf, false);
 
     psTimerStart ("psphot.apresid");
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    assert (recipe);
+
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
+    // nThreads = 0; // XXX until testing is complete, do not thread this function
 
     bool measureAptrend = psMetadataLookupBool (&status, recipe, "MEASURE.APTREND");
@@ -70 +81 @@
     pmSourceMagnitudesInit (recipe);
 
+    // threaded measurement of the source magnitudes
+    // 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_APRESID_MAGS");
+
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, psf);
+            PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_U8); // XXX change this to use abstract mask type info
+
+            PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nskip
+            PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nfail
+
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree (job);
+                return false;
+            }
+            psFree(job);
+
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            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];
+                Nskip += scalar->data.S32;
+                scalar = job->args->data[5];
+                Nfail += scalar->data.S32;
+            }
+            psFree(job);
+        }
+    }
+
+    psFree (cellGroups);
+
+    // gather the stats to assess the aperture residuals
     psVector *mask    = psVectorAllocEmpty (300, PS_TYPE_U8);
     psVector *mag     = psVectorAllocEmpty (300, PS_TYPE_F32);
@@ -78 +147 @@
     Npsf = 0;
 
-    // select all good PM_SOURCE_TYPE_STAR entries
     for (int i = 0; i < sources->n; i++) {
         source = sources->data[i];
@@ -89 +157 @@
         if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
 
-        // get growth-corrected, apTrend-uncorrected magnitudes in scaled apertures
-        // will fail if below S/N threshold or model is missing
-        if (!pmSourceMagnitudes (source, psf, photMode, maskVal)) {
-            Nskip ++;
-            psTrace ("psphot", 3, "skip : bad source mag");
-            continue;
-        }
-
         if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
-            Nfail ++;
-            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
             continue;
         }
@@ -206 +264 @@
 
     psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);
-    psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for: %f sec\n", psTimerMark ("psphot.apresid"));
+    psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d sources: %f sec\n", Npsf, psTimerMark ("psphot.apresid"));
 
     psFree (mag);
@@ -392 +450 @@
 }
 
+bool psphotApResidMags_Threaded (psThreadJob *job) {
+
+    int Nskip = 0;
+    int Nfail = 0;
+
+    psScalar *scalar = NULL;
+
+    psArray *sources                = job->args->data[0];
+    pmPSF *psf                      = job->args->data[1];
+    pmSourcePhotometryMode photMode = PS_SCALAR_VALUE(job->args->data[2],S32);
+    psMaskType maskVal              = PS_SCALAR_VALUE(job->args->data[3],U8);
+
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = (pmSource *) sources->data[i];
+
+        if (source->type != PM_SOURCE_TYPE_STAR) SKIPSTAR ("NOT STAR");
+        if (source->mode &  PM_SOURCE_MODE_SATSTAR) SKIPSTAR ("SATSTAR");
+        if (source->mode &  PM_SOURCE_MODE_BLEND) SKIPSTAR ("BLEND");
+        if (source->mode &  PM_SOURCE_MODE_FAIL) SKIPSTAR ("FAIL STAR");
+        if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
+
+        if (!pmSourceMagnitudes (source, psf, photMode, maskVal)) {
+            Nskip ++;
+            psTrace ("psphot", 3, "skip : bad source mag");
+            continue;
+        }
+
+        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
+            Nfail ++;
+            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
+            continue;
+        }
+    }
+
+    // change the value of a scalar on the array (wrap this and put it in psArray.h)
+    scalar = job->args->data[4];
+    scalar->data.S32 = Nskip;
+
+    scalar = job->args->data[5];
+    scalar->data.S32 = Nfail;
+
+    return true;
+}

trunk/psphot/src/psphotArguments.c

@@ -38 +38 @@
 
         // create the thread pool with number of desired threads, supplying our thread launcher function
-        // XXX need to determine the number of threads from the config data
         psThreadPoolInit (nThreads);
     }

trunk/psphot/src/psphotBlendFit.c

@@ -2 +2 @@
 
 // XXX I don't like this name
-bool psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) {
+bool psphotBlendFit (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) {
 
     int Nfit = 0;
@@ -12 +12 @@
     psTimerStart ("psphot.fit.nonlinear");
 
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    assert (recipe);
+
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
+    // nThreads = 0; // XXX until testing is complete, do not thread this function
+
+    psphotInitLimitsPSF (recipe, readout);
+    psphotInitLimitsEXT (recipe);
+    psphotInitRadiusPSF (recipe, psf->type);
+
+    // starts the timer, sets up the array of fitSets 
+    psphotFitInit (nThreads);
+
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, pmSourceSortBySN);
+
+    // 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);
+
+    fprintf (stderr, "starting with %ld sources\n", sources->n);
+
+    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_BLEND_FIT");
+            psArray *newSources = psArrayAllocEmpty(16);
+
+            psArrayAdd(job->args, 1, readout);
+            psArrayAdd(job->args, 1, recipe);
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, psf);
+            psArrayAdd(job->args, 1, newSources); // return for new sources
+            psFree (newSources);
+
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf
+            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 Nfail
+
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree (job);
+                return NULL;
+            }
+            psFree(job);
+
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false)) {
+            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 = NULL;
+                scalar = job->args->data[5];
+                Nfit += scalar->data.S32;
+                scalar = job->args->data[6];
+                Npsf += scalar->data.S32;
+                scalar = job->args->data[7];
+                Next += scalar->data.S32;
+                scalar = job->args->data[8];
+                Nfail += scalar->data.S32;
+
+                // add these back onto sources
+                psArray *newSources = job->args->data[4];
+                for (int j = 0; j < newSources->n; j++) {
+                    psArrayAdd (sources, 16, newSources->data[j]);
+                }
+            }
+            psFree(job);
+        }
+    }
+    psFree (cellGroups);
+
+    if (psTraceGetLevel("psphot") >= 6) {
+      psphotSaveImage (NULL, readout->image,  "image.v2.fits");
+    }
+
+    psLogMsg ("psphot.psphotBlendFit", PS_LOG_INFO, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot.fit.nonlinear"), Nfit, Npsf, Next, Nfail, sources->n - Nfit);
+
+    psphotVisualShowResidualImage (readout);
+    psphotVisualShowFlags (sources);
+
+    return true;
+}
+
+bool psphotBlendFit_Threaded (psThreadJob *job) {
+
+    bool status = false;
+    int Nfit = 0;
+    int Npsf = 0;
+    int Next = 0;
+    int Nfail = 0;
+    psScalar *scalar = NULL;
+
+    pmReadout *readout  = job->args->data[0];
+    psMetadata *recipe  = job->args->data[1];
+    psArray *sources    = job->args->data[2];
+    pmPSF *psf          = job->args->data[3];
+    psArray *newSources = job->args->data[4];
+
     // bit-masks to test for good/bad pixels
     psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT");
@@ -23 +143 @@
     maskVal |= markVal;
 
-    // source analysis is done in S/N order (brightest first)
-    sources = psArraySort (sources, pmSourceSortBySN);
-
     // S/N limit to perform full non-linear fits
     float FIT_SN_LIM = psMetadataLookupF32 (&status, recipe, "FULL_FIT_SN_LIM");
-
-    psphotInitLimitsPSF (recipe, readout);
-    psphotInitLimitsEXT (recipe);
-    psphotInitRadiusPSF (recipe, psf->type);
-
-    psphotFitInit ();
 
     // option to limit analysis to a specific region
@@ -41 +152 @@
 
     for (int i = 0; i < sources->n; i++) {
-        // if (i%100 == 0) psphotFitSummary ();
-
         pmSource *source = sources->data[i];
 
@@ -57 +166 @@
         if (source->peak->SN < FIT_SN_LIM) continue;
 
-        // XXX this should use peak?
+        // exclude sources outside optional analysis region
         if (source->peak->xf < AnalysisRegion.x0) continue;
         if (source->peak->yf < AnalysisRegion.y0) continue;
@@ -66 +175 @@
         if (source->modelPSF == NULL) continue;
 
-        // skip sources which are insignificant flux?
+        // skip sources which are insignificant flux? 
+        // XXX this is somewhat ad-hoc
         if (source->modelPSF->params->data.F32[1] < 0.1) {
             psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n",
@@ -81 +191 @@
         Nfit ++;
 
-        // XXX TEST
-        // if ((fabs(source->peak->x - 1202) < 2) && (fabs(source->peak->y - 1065) < 2)) {
-        // psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6);
-        // }
-
-        // try fitting PSFs, then try extended sources
-        // these functions subtract the resulting fitted source (XXX and update the modelFlux?)
-        // XXX re-consider conditions under which the source has EXT fit:
-        // I should try EXT if the source size measurement says it is large
+        // try fitting PSFs or extended sources depending on source->mode
+        // these functions subtract the resulting fitted source
         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
-            if (psphotFitBlob (readout, source, sources, psf, maskVal, markVal)) {
+            if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
                 source->type = PM_SOURCE_TYPE_EXTENDED;
                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
@@ -115 +218 @@
     }
 
-    if (psTraceGetLevel("psphot") >= 6) {
-      psphotSaveImage (NULL, readout->image,  "image.v2.fits");
-    }
-
-    psLogMsg ("psphot.psphotBlendFit", PS_LOG_INFO, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot.fit.nonlinear"), Nfit, Npsf, Next, Nfail, sources->n - Nfit);
-
-    psphotVisualShowResidualImage (readout);
-    psphotVisualShowFlags (sources);
+    // 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 = Nfit;
+
+    scalar = job->args->data[6];
+    scalar->data.S32 = Npsf;
+
+    scalar = job->args->data[7];
+    scalar->data.S32 = Next;
+
+    scalar = job->args->data[8];
+    scalar->data.S32 = Nfail;
 
     return true;

trunk/psphot/src/psphotCleanup.c

@@ -12 +12 @@
     pmConceptsDone ();
     pmConfigDone ();
+    pmSourceFitSetDone ();
     // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psphot");
     fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psphot");

trunk/psphot/src/psphotGuessModels.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-bool psphotChooseCellSizes (int *Cx, int *Cy, pmReadout *readout, int nThreads);
-psphotGuessModelForRegionArgs *psphotGuessModelForRegionArgsAlloc();
-bool psphotGuessModelForRegion (psphotGuessModelForRegionArgs *args);
+// XXX : the threading here is not great.  this may be due to blocks between elements, but
+// the selection of the objects in a cell is not optimal.  To fix:
+// 1) define the boundaries of the cells up front 
+// 2) loop over the sources once and associate them with their cell
+// 3) define the threaded function to work with sources for a given cell
 
 // A guess for when the moments aren't available
@@ -20 +22 @@
 }
 
-# define NFILL 4
-
 // construct an initial PSF model for each object
 bool psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) {
@@ -33 +33 @@
     assert (recipe);
 
-    // int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
-    // if (!status) {
-    // nThreads = 0;
-    // }
-    int nThreads = 0;
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
+    // nThreads = 0; // XXX until testing is complete, do not thread this function
 
     // bit-masks to test for good/bad pixels
@@ -53 +54 @@
     psphotInitRadiusPSF (recipe, psf->type);
 
-    // the strategy here is to divide the image into 2x2 blocks of cells and cycle through
-    // the four discontiguous sets of cells, threading all within a set and blocking between
-    // sets 
-
-    // we divide the image region into 2*2 blocks of size Nx*Ny, the image will have 
-    // Cx*Cy blocks so that (2Nx)Cx = numCols, (2Ny)Cy = numRows.  We want to choose Cx and
-    // Cy so that (2Nx)Cx * (2Ny)Cy = 4 * NFILL * nThreads -- each of the four sets of cells
-    // has enough cells to allow NFILL cells for each thread (to better distribute heavy and
-    // light load cells
-    
-    // choose Cx, Cy:
+    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
     int Cx = 1, Cy = 1;
     psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
 
-    // the array runs from readout->image->col0 to readout->image->col0 + readout->image->numCols 
-    int Xo = readout->image->col0;
-    int Yo = readout->image->row0;
-    int Nx = readout->image->numCols / (2*Cx);
-    int Ny = readout->image->numRows / (2*Cy);
-
-    // we can thread all of the cells within a set, but need to block between sets
-    for (int jy = 0; jy < 2; jy++) {
-        for (int jx = 0; jx < 2; jx++) {
-
-            // generate jobs for all of the Cx*Cy cells in this set
-            for (int ix = 0; ix < Cx; ix++) {
-                for (int iy = 0; iy < Cy; iy++) {
-                
-                    int x0 = (2*ix + jx)*Nx + Xo;
-                    int x1 = x0 + Nx;
-
-                    int y0 = (2*iy + jy)*Ny + Yo;
-                    int y1 = y0 + Ny;
-
-                    if (readout->image->numCols + Xo - x1 < Nx) {
-                        x1 = readout->image->numCols + Xo;
-                    }
-                    if (readout->image->numRows + Yo - y1 < Ny) {
-                        y1 = readout->image->numRows + Yo;
-                    }
-
-                    // fprintf (stderr, "launch %d,%d - %d,%d\n", x0, y0, x1, y1);
-                    
-                    psRegion *cellRegion = psRegionAlloc (x0, x1, y0, y1);
-                    *cellRegion = psRegionForImage (readout->image, *cellRegion);
-                    // if we got the math above right, this should be a NOP
-                                                                             
-                    psphotGuessModelForRegionArgs *args = psphotGuessModelForRegionArgsAlloc();
-                    args->readout = psMemIncrRefCounter(readout);
-                    args->sources = psMemIncrRefCounter(sources);
-                    args->psf     = psMemIncrRefCounter(psf);
-                    args->region  = psMemIncrRefCounter(cellRegion);
-
-                    args->maskVal = maskVal;
-                    args->markVal = markVal;
-
-                    // allocate a job -- if threads are not defined, this just runs the job
-                    psThreadJob *job = psThreadJobAlloc ("PSPHOT_GUESS_MODEL");
-                    psArrayAdd(job->args, 1, args);
-                    if (!psThreadJobAddPending(job)) {
-                        psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                        return false;
-                    }
-                    psFree(cellRegion);
-                    psFree(job);
-                    psFree(args);
-                }
-            }
-
-            // wait for the threads to finish and manage results
-            // wait here for the threaded jobs to finish
-            // fprintf (stderr, "wait for threads (%d, %d)\n", jx, jy);
-            if (!psThreadPoolWait (false)) {
+    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_GUESS_MODEL");
+            psArrayAdd(job->args, 1, readout);
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, psf);
+
+            // XXX change these to use abstract mask type info
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_U8);
+            PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_U8);
+
+            if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree (job);
                 return false;
             }
-
-            // we have only supplied one type of job, so we can assume the types here
-            psThreadJob *job = NULL;
-            while ((job = psThreadJobGetDone()) != NULL) {
-                // we have no returned data from this operation
-                if (job->args->n < 1) {
-                    fprintf (stderr, "error with job\n");
-                }
-                psFree(job);
-            }
-        }
-    }
-    
+            psFree(job);
+        }
+
+        // wait for the threads to finish and manage results
+        // wait here for the threaded jobs to finish
+        // fprintf (stderr, "wait for threads (%d, %d)\n", jx, jy);
+        if (!psThreadPoolWait (false)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            return false;
+        }
+
+        // we have only supplied one type of job, so we can assume the types here
+        psThreadJob *job = NULL;
+        while ((job = psThreadJobGetDone()) != NULL) {
+            // we have no returned data from this operation
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            }
+            psFree(job);
+        }
+    }
+
+    // XXX check that all sources that should have models
+    int nMiss = 0;
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+            source->mode &= ~PM_SOURCE_MODE_EXT_LIMIT;
+            continue;
+        }
+
+        nMiss ++;
+    }
+    psLogMsg ("psphot.models", 4, "failed to build models for %d objects\n", nMiss);
+
+    psFree (cellGroups);
+
     psLogMsg ("psphot.models", 4, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.models"));
     return true;
 }
 
-static void psphotGuessModelForRegionArgsFree(psphotGuessModelForRegionArgs *args)
-{
-    psFree(args->readout);
-    psFree(args->sources);
-    psFree(args->psf);
-    psFree(args->region);
-    return;
-}
-
-psphotGuessModelForRegionArgs *psphotGuessModelForRegionArgsAlloc()
-{
-    psphotGuessModelForRegionArgs *args = psAlloc(sizeof(psphotGuessModelForRegionArgs));
-    psMemSetDeallocator(args, (psFreeFunc)psphotGuessModelForRegionArgsFree);
-
-    args->readout = NULL;
-    args->sources = NULL;
-    args->psf = NULL;
-    args->region = NULL;
-
-    args->maskVal = 0;
-    args->markVal = 0;
-    return args;
-}
-
 // construct models only for sources in the specified region
-bool psphotGuessModelForRegion (psphotGuessModelForRegionArgs *args) {
-
-    pmReadout *readout = args->readout;
-    psArray *sources = args->sources;
-    pmPSF *psf = args->psf;
-    psRegion *region = args->region;
-    psMaskType maskVal = args->maskVal;
-    psMaskType markVal = args->markVal;
+bool psphotGuessModel_Threaded (psThreadJob *job) {
+
+    pmReadout *readout = job->args->data[0];
+    psArray *sources   = job->args->data[1];
+    pmPSF *psf         = job->args->data[2];
+    
+    psMaskType maskVal = PS_SCALAR_VALUE(job->args->data[3],U8);
+    psMaskType markVal = PS_SCALAR_VALUE(job->args->data[4],U8);
 
     int nSrc = 0;
@@ -184 +137 @@
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = sources->data[i];
+
+        // XXXX this is just for a test: use this to mark sources for which the model is measured
+        // check later that all are used.
+        source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
 
         // skip non-astronomical objects (very likely defects)
@@ -189 +146 @@
         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
         if (!source->peak) continue;
-
-        if (source->peak->xf <  region->x0) continue;
-        if (source->peak->yf <  region->y0) continue;
-        if (source->peak->xf >= region->x1) continue;
-        if (source->peak->yf >= region->y1) continue;
 
         nSrc ++;
@@ -256 +208 @@
     }
 
-    // fprintf (stderr, "%d for region %lf,%lf - %lf,%lf\n", nSrc, region->x0, region->y0, region->x1, region->y1);
     return true;
 }
 
-bool psphotChooseCellSizes (int *Cx, int *Cy, pmReadout *readout, int nThreads) {
-
-    int nCells = nThreads * NFILL; // number of cells in a single set
-    int C = sqrt(nCells) + 0.5;
-    
-    // we need to assign Cx and Cy based on the dimensionality of the image
-    // crude way to find most evenly balanced factors of nCells:
-    for (int i = C; i >= 1; i--) {
-        int C1 = nCells / C;
-        int C2 = nCells / C1;
-        if (C1*C2 != nCells) continue;
-
-        if (readout->image->numRows > readout->image->numCols) {
-            *Cx = PS_MAX (C1, C2);
-            *Cy = PS_MIN (C1, C2);
-        } else {
-            *Cx = PS_MAX (C1, C2);
-            *Cy = PS_MIN (C1, C2);
-        }
-        return true;
-    }
-    *Cx = 1;
-    *Cy = 1; 
-
-    return true;
-}

trunk/psphot/src/psphotMagnitudes.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-bool psphotMagnitudes(pmConfig *config, const pmFPAview *view, psArray *sources, psMetadata *recipe, pmPSF *psf) {
+bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf) {
 
     bool status = false;
@@ -7 +7 @@
 
     psTimerStart ("psphot.mags");
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    assert (recipe);
+
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
+    // nThreads = 0; // XXX until testing is complete, do not thread this function
 
     pmReadout *backModel = psphotSelectBackground (config, view);
@@ -40 +51 @@
     if (INTERPOLATE_AP) photMode |= PM_SOURCE_PHOT_INTERP;
 
+    // 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_MAGNITUDES");
+
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, psf);
+            psArrayAdd(job->args, 1, binning);
+            psArrayAdd(job->args, 1, backModel);
+            psArrayAdd(job->args, 1, backStdev);
+
+            PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_U8); // XXX change this to use abstract mask type info
+            PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for nAp
+
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree (job);
+                return false;
+            }
+            psFree(job);
+
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            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 = job->args->data[7];
+                Nap += scalar->data.S32;
+            }
+            psFree(job);
+        }
+    }
+
+    psFree (cellGroups);
+
+    psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "measure magnitudes : %f sec for %ld objects (%d with apertures)\n", psTimerMark ("psphot.mags"), sources->n, Nap);
+    return true;
+}
+
+bool psphotMagnitudes_Threaded (psThreadJob *job) {
+
+    bool status;
+    int Nap = 0;
+
+    psArray *sources                = job->args->data[0];
+    pmPSF *psf                      = job->args->data[1];
+    psImageBinning *binning         = job->args->data[2];
+    pmReadout *backModel            = job->args->data[3];
+    pmReadout *backStdev            = job->args->data[4];
+    pmSourcePhotometryMode photMode = PS_SCALAR_VALUE(job->args->data[5],S32);
+    psMaskType maskVal              = PS_SCALAR_VALUE(job->args->data[6],U8);
+
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
@@ -47 +130 @@
         source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, backModel->image, binning);
         if (isnan(source->sky) && false) {
-          psError(PSPHOT_ERR_SKY, false, "Setting pmSource.sky");
-          psErrorStackPrint(NULL, " ");
-          psErrorClear();
+            psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.sky");
         }
 
         source->skyErr = psImageUnbinPixel(source->peak->x, source->peak->y, backStdev->image, binning);
         if (isnan(source->skyErr) && false) {
-          psError(PSPHOT_ERR_SKY, false, "Setting pmSource.sky");
-          psErrorStackPrint(NULL, " ");
-          psErrorClear();
+            psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.skyErr");
         }
     }
 
-    psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "measure magnitudes : %f sec for %ld objects (%d with apertures)\n", psTimerMark ("psphot.mags"), sources->n, Nap);
+    // change the value of a scalar on the array (wrap this and put it in psArray.h)
+    psScalar *scalar = job->args->data[7];
+    scalar->data.S32 = Nap;
+
     return true;
 }
-
-// XXX add in a measurement of the bright and faint completeness values
-// XXX push these into the recipe

trunk/psphot/src/psphotReadout.c

@@ -79 +79 @@
 
     // construct sources and measure basic stats
-    psArray *sources = psphotSourceStats (readout, recipe, detections);
+    psArray *sources = psphotSourceStats (config, readout, detections);
     if (!sources) return false;
     if (!strcasecmp (breakPt, "PEAKS")) {
@@ -178 +178 @@
 
     // non-linear PSF and EXT fit to brighter sources
-    psphotBlendFit (readout, sources, recipe, psf);
+    psphotBlendFit (config, readout, sources, psf);
 
     // replace all sources
@@ -207 +207 @@
 
     // define new sources based on only the new peaks
-    psArray *newSources = psphotSourceStats (readout, recipe, detections);
+    psArray *newSources = psphotSourceStats (config, readout, detections);
 
     // set source type
@@ -243 +243 @@
 
     // measure aperture photometry corrections
-    if (!psphotApResid (readout, sources, recipe, psf)) {
+    if (!psphotApResid (config, readout, sources, psf)) {
         psLogMsg ("psphot", 3, "failed on psphotApResid");
         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
@@ -249 +249 @@
 
     // calculate source magnitudes
-    psphotMagnitudes(config, view, sources, recipe, psf);
+    psphotMagnitudes(config, readout, view, sources, psf);
 
     // replace failed sources?

trunk/psphot/src/psphotReadoutFindPSF.c

@@ -37 +37 @@
 
     // construct sources and measure basic stats (moments, local sky)
-    psArray *sources = psphotSourceStats(readout, recipe, detections);
+    psArray *sources = psphotSourceStats(config, readout, detections);
     if (!sources) return false;
 

trunk/psphot/src/psphotSetThreads.c

@@ -1 +1 @@
 # include "psphot.h"
-
-// each thread runs this function, starting a new job when it finished with an old one
-// it is called with a (void *) pointer to its own thread pointer
-bool psphotThread_psphotGuessModel(psThreadJob *job)
-{
-    psphotGuessModelForRegionArgs *args = job->args->data[0];
-    bool status = psphotGuessModelForRegion (args);
-    return status;
-}
 
 bool psphotSetThreads () {
@@ -14 +5 @@
     psThreadTask *task = NULL;
 
-    task = psThreadTaskAlloc("PSPHOT_GUESS_MODEL", 1);
-    task->function = &psphotThread_psphotGuessModel;
+    task = psThreadTaskAlloc("PSPHOT_GUESS_MODEL", 5);
+    task->function = &psphotGuessModel_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    task = psThreadTaskAlloc("PSPHOT_MAGNITUDES", 8);
+    task->function = &psphotMagnitudes_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    task = psThreadTaskAlloc("PSPHOT_APRESID_MAGS", 6);
+    task->function = &psphotApResidMags_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    task = psThreadTaskAlloc("PSPHOT_SOURCE_STATS", 4);
+    task->function = &psphotSourceStats_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    task = psThreadTaskAlloc("PSPHOT_BLEND_FIT", 9);
+    task->function = &psphotBlendFit_Threaded;
     psThreadTaskAdd(task);
     psFree(task);

trunk/psphot/src/psphotSourceFits.c

@@ -9 +9 @@
 static int NfitEXT = 0;
 
-bool psphotFitInit () {
+bool psphotFitInit (int nThreads) {
     psTimerStart ("psphot.fits");
+    pmSourceFitSetInit (nThreads);
     return true;
 }
@@ -207 +208 @@
 }
 
-bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, psMaskType maskVal, psMaskType markVal) {
+bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, psMaskType maskVal, psMaskType markVal) {
 
     bool okEXT, okDBL;
@@ -310 +311 @@
     psTrace ("psphot", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[PM_PAR_XPOS], ONE->params->data.F32[PM_PAR_YPOS]);
 
-    psArrayAdd (sources, 100, newSrc);
+    psArrayAdd (newSources, 100, newSrc);
     psFree (newSrc);
     psFree (DBL);

trunk/psphot/src/psphotSourceStats.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, pmDetections *detections)
-{
-    bool     status  = false;
+psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections) {
+
+    bool status = false;
     psArray *sources = NULL;
-    float BIG_RADIUS;
 
     psTimerStart ("psphot.stats");
 
-    // bit-masks to test for good/bad pixels
-    psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT");
-    assert (maskVal);
-
-    // bit-mask to mark pixels not used in analysis
-    psMaskType markVal = psMetadataLookupU8(&status, recipe, "MARK.PSPHOT");
-    assert (markVal);
-
-    // maskVal is used to test for rejected pixels, and must include markVal
-    maskVal |= markVal;
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    assert (recipe);
+
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
+    // nThreads = 0; // XXX until testing is complete, do not thread this function
 
     // determine properties (sky, moments) of initial sources
-    float INNER    = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
-    if (!status) return NULL;
     float OUTER    = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
-    if (!status) return NULL;
-    float RADIUS   = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
-    if (!status) return NULL;
-    float MIN_SN   = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
     if (!status) return NULL;
     char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
@@ -38 +31 @@
     }
 
+    // generate the array of sources, define the associated pixel
     sources = psArrayAllocEmpty (peaks->n);
 
+    for (int i = 0; i < peaks->n; i++) {
+
+        pmPeak *peak = peaks->data[i];
+        if (peak->assigned) continue;
+
+        // create a new source
+        pmSource *source = pmSourceAlloc();
+
+        // add the peak
+        source->peak = psMemIncrRefCounter(peak);
+
+        // allocate space for moments
+        source->moments = pmMomentsAlloc();
+
+        // allocate image, weight, mask arrays for each peak (square of radius OUTER)
+        pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
+
+        peak->assigned = true;
+        psArrayAdd (sources, 100, source);
+        psFree (source);
+    }
+
+    if (!strcasecmp (breakPt, "PEAKS")) {
+        psLogMsg ("psphot", PS_LOG_INFO, "%ld sources : %f sec\n", sources->n, psTimerMark ("psphot.stats"));
+        psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n");
+        psphotVisualShowMoments (sources);
+        return sources;
+    }
+
+    // threaded measurement of the source magnitudes
     int Nfail = 0;
     int Nmoments = 0;
-    for (int i = 0; i < peaks->n; i++) {
-
-        pmPeak *peak = peaks->data[i];
-        if (peak->assigned) continue;
-
-        // create a new source, add peak
-        pmSource *source = pmSourceAlloc();
-        source->peak = psMemIncrRefCounter(peak);
-
-        // allocate image, weight, mask arrays for each peak (square of radius OUTER)
-        pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
-        if (!strcasecmp (breakPt, "PEAKS")) {
-            peak->assigned = true;
-            psArrayAdd (sources, 100, source);
-            psFree (source);
-            continue;
-        }
+
+    // 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_SOURCE_STATS");
+
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, recipe);
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
+
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree (job);
+                return NULL;
+            }
+            psFree(job);
+
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false)) {
+            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 = NULL;
+                scalar = job->args->data[2];
+                Nmoments += scalar->data.S32;
+                scalar = job->args->data[3];
+                Nfail += scalar->data.S32;
+            }
+            psFree(job);
+        }
+    }
+
+    psFree (cellGroups);
+
+    psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d failed: %f sec\n", sources->n, Nmoments, Nfail, psTimerMark ("psphot.stats"));
+
+    psphotVisualShowMoments (sources);
+
+    return (sources);
+}
+
+bool psphotSourceStats_Threaded (psThreadJob *job) {
+
+    bool status = false;
+    float BIG_RADIUS;
+    psScalar *scalar = NULL;
+
+    psArray *sources                = job->args->data[0];
+    psMetadata *recipe              = job->args->data[1];
+
+    float INNER    = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
+    if (!status) return false;
+    float MIN_SN   = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
+    if (!status) return false;
+    float RADIUS   = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
+    if (!status) return false;
+
+    // bit-masks to test for good/bad pixels
+    psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT");
+    assert (maskVal);
+
+    // bit-mask to mark pixels not used in analysis
+    psMaskType markVal = psMetadataLookupU8(&status, recipe, "MARK.PSPHOT");
+    assert (markVal);
+
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
+
+    // threaded measurement of the sources moments
+    int Nfail = 0;
+    int Nmoments = 0;
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
 
         // skip faint sources for moments measurement
         if (source->peak->SN < MIN_SN) {
-            peak->assigned = true;
-            psArrayAdd (sources, 100, source);
-            psFree (source);
             continue;
         }
@@ -72 +169 @@
         status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
         if (!status) {
-          psFree (source);
-          Nfail ++;
-          psErrorClear();
-          continue;
+            psErrorClear(); // XXX re-consider the errors raised here
+            Nfail ++;
+            continue;
         }
 
@@ -82 +178 @@
         status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
         if (!status) {
-          psFree (source);
-          Nfail ++;
-          psErrorClear();
-          continue;
+            Nfail ++;
+            psErrorClear();
+            continue;
         }
 
@@ -91 +186 @@
         status = pmSourceMoments (source, RADIUS);
         if (status) {
-            // add to the source array
-            peak->assigned = true;
-            psArrayAdd (sources, 100, source);
-            psFree (source);
             Nmoments ++;
             continue;
@@ -105 +196 @@
         status = pmSourceMoments (source, BIG_RADIUS);
         if (status) {
-            // add to the source array
-            peak->assigned = true;
-            psArrayAdd (sources, 100, source);
-            psFree (source);
             Nmoments ++;
             continue;
         }
 
-        psFree (source);
         Nfail ++;
         psErrorClear();
@@ -119 +205 @@
     }
 
-    psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d failed: %f sec\n", sources->n, Nmoments, Nfail, psTimerMark ("psphot.stats"));
-
-    psphotVisualShowMoments (sources);
-
-    return (sources);
+    // change the value of a scalar on the array (wrap this and put it in psArray.h)
+    scalar = job->args->data[2];
+    scalar->data.S32 = Nmoments;
+
+    scalar = job->args->data[3];
+    scalar->data.S32 = Nfail;
+    
+    return true;
 }
-
-// XXX EAM : filter out bad peaks (eg, no valid pixels available for sky)