Index: trunk/psphot/src/psphotBlendFit.c
===================================================================
--- trunk/psphot/src/psphotBlendFit.c	(revision 21250)
+++ trunk/psphot/src/psphotBlendFit.c	(revision 21359)
@@ -1,3 +1,5 @@
 # include "psphotInternal.h"
+
+bool psphotBlendFit_Unthreaded (int *nfit, int *npsf, int *next, int *nfail, pmReadout *readout, psMetadata *recipe, psArray *sources, pmPSF *psf, psArray *newSources);
 
 // XXX I don't like this name
@@ -21,5 +23,4 @@
 	nThreads = 0;
     }
-    // nThreads = 0; // XXX until testing is complete, do not thread this function
 
     // source fitting parameters for extended source fits
@@ -63,58 +64,82 @@
 	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)) {
+	    if (nThreads) {
+		// 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);
+	    } else {
+		int nfit = 0;
+		int npsf = 0;
+		int next = 0;
+		int nfail = 0;
+		psArray *newSources = psArrayAllocEmpty(16);
+
+		if (!psphotBlendFit_Unthreaded (&nfit, &npsf, &next, &nfail, readout, recipe, cells->data[j], psf, newSources)) {
+		    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+		    return NULL;
+		}
+		Nfit += nfit;
+		Npsf += npsf;
+		Next += next;
+		Nfail += nfail;
+
+		// add these back onto sources
+		for (int k = 0; k < newSources->n; k++) {
+		    psArrayAdd (sources, 16, newSources->data[k]);
+		}
+		psFree (newSources);
+	    }
+	}
+
+	if (nThreads) {
+	    // wait for the threads to finish and manage results
+	    if (!psThreadPoolWait (false)) {
 		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]);
+
+	    // 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(job);
+	    }
 	}
     }
@@ -249,2 +274,105 @@
     return true;
 }
+
+bool psphotBlendFit_Unthreaded (int *nfit, int *npsf, int *next, int *nfail, pmReadout *readout, psMetadata *recipe, psArray *sources, pmPSF *psf, psArray *newSources) {
+
+    bool status = false;
+    int Nfit = 0;
+    int Npsf = 0;
+    int Next = 0;
+    int Nfail = 0;
+
+    // bit-masks to test for good/bad pixels
+    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
+    assert (maskVal);
+
+    // bit-mask to mark pixels not used in analysis
+    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
+    assert (markVal);
+
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
+
+    // S/N limit to perform full non-linear fits
+    float FIT_SN_LIM = psMetadataLookupF32 (&status, recipe, "FULL_FIT_SN_LIM");
+
+    // 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");
+
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+
+        // skip non-astronomical objects (very likely defects)
+        if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
+        if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) continue;
+        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+
+        // skip DBL second sources (ie, added by psphotFitBlob)
+        if (source->mode &  PM_SOURCE_MODE_PAIR) continue;
+
+        // limit selection to some SN limit
+        if (source->peak->SN < FIT_SN_LIM) continue;
+
+	// exclude sources outside optional analysis region
+        if (source->peak->xf < AnalysisRegion.x0) continue;
+        if (source->peak->yf < AnalysisRegion.y0) continue;
+        if (source->peak->xf > AnalysisRegion.x1) continue;
+        if (source->peak->yf > AnalysisRegion.y1) continue;
+
+        // if model is NULL, we don't have a starting guess
+        if (source->modelPSF == NULL) continue;
+
+        // 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",
+                     source->modelPSF->params->data.F32[1],
+                     source->modelPSF->params->data.F32[2],
+                     source->modelPSF->params->data.F32[3]);
+            continue;
+        }
+
+        // replace object in image
+        if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
+            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+        }
+        Nfit ++;
+
+        // 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, 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);
+		Next ++;
+		continue;
+	    }
+	} else {
+	    if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
+		source->type = PM_SOURCE_TYPE_STAR;
+		psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
+		Npsf ++;
+		continue;
+	    }
+	}
+
+        psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf);
+        Nfail ++;
+
+        // re-subtract the object, leave local sky
+        pmSourceCacheModel (source, maskVal);
+        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+        source->mode |= PM_SOURCE_MODE_SUBTRACTED;
+    }
+
+    // change the value of a scalar on the array (wrap this and put it in psArray.h)
+    *nfit  = Nfit;
+    *npsf  = Npsf;
+    *next  = Next;
+    *nfail = Nfail;
+
+    return true;
+}