Index: trunk/psphot/src/psphotGuessModels.c
===================================================================
--- trunk/psphot/src/psphotGuessModels.c	(revision 20453)
+++ trunk/psphot/src/psphotGuessModels.c	(revision 21166)
@@ -1,7 +1,9 @@
 # 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,6 +22,4 @@
 }
 
-# define NFILL 4
-
 // construct an initial PSF model for each object
 bool psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) {
@@ -33,9 +33,10 @@
     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,130 +54,82 @@
     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,4 +137,8 @@
     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,9 +146,4 @@
 	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,32 +208,5 @@
     }
 
-    // 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;
-}