Index: trunk/psphot/src/psphotSourceStats.c
===================================================================
--- trunk/psphot/src/psphotSourceStats.c	(revision 20453)
+++ trunk/psphot/src/psphotSourceStats.c	(revision 21166)
@@ -1,31 +1,24 @@
 # 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 +31,135 @@
     }
 
+    // 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,8 +169,7 @@
         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,8 +178,7 @@
         status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
         if (!status) {
-          psFree (source);
-          Nfail ++;
-          psErrorClear();
-          continue;
+	    Nfail ++;
+	    psErrorClear();
+	    continue;
         }
 
@@ -91,8 +186,4 @@
         status = pmSourceMoments (source, RADIUS);
         if (status) {
-            // add to the source array
-            peak->assigned = true;
-            psArrayAdd (sources, 100, source);
-            psFree (source);
             Nmoments ++;
             continue;
@@ -105,13 +196,8 @@
         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,10 +205,11 @@
     }
 
-    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)