Index: trunk/psphot/src/psphotExtendedSourceAnalysis.c
===================================================================
--- trunk/psphot/src/psphotExtendedSourceAnalysis.c	(revision 31154)
+++ trunk/psphot/src/psphotExtendedSourceAnalysis.c	(revision 32348)
@@ -2,5 +2,4 @@
 
 // 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
@@ -8,4 +7,7 @@
 {
     bool status = true;
+
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Extended Source Analysis (Petrosians) ---");
 
     // select the appropriate recipe information
@@ -59,7 +61,9 @@
     }
 
-    // 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,4 +74,120 @@
     }
 
+    // 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_WARN, "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,11 +198,7 @@
     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
@@ -90,4 +206,9 @@
 
 	pmSource *source = sources->data[i];
+
+	// if we have checked the source validity on the basis of the object set, then 
+	// we either skip these tests below or we skip the source completely
+	if (source->tmpFlags & PM_SOURCE_TMPF_STACK_SKIP) continue;
+	if (source->tmpFlags & PM_SOURCE_TMPF_STACK_KEEP) goto keepSource;
 
 	// skip PSF-like and non-astronomical objects
@@ -104,12 +225,21 @@
 
 	// limit selection to some SN limit
-	assert (source->peak); // how can a source not have a peak?
-	if (sqrt(source->peak->detValue) < SN_LIM) continue;
-
-	// 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;
+	// assert (source->peak); // how can a source not have a peak?
+	// limit selection to some SN limit
+	bool skipSource = false;
+	if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
+	    skipSource = (source->moments->KronFlux < SN_LIM * source->moments->KronFluxErr);
+	} else {
+	    skipSource = (sqrt(source->peak->detValue) < SN_LIM);
+	}
+	if (skipSource) continue;
+
+	// limit selection by analysis region (this automatically apply
+	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;
+
+    keepSource:
 
 	// replace object in image
@@ -159,15 +289,15 @@
     }
 
-    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;