Index: trunk/psphot/src/Makefile.am =================================================================== --- trunk/psphot/src/Makefile.am (revision 20356) +++ trunk/psphot/src/Makefile.am (revision 20411) @@ -18,4 +18,5 @@ libpsphot_la_SOURCES = \ + psphotSetThreads.c \ psphotErrorCodes.c \ psphotImageQuality.c \ Index: trunk/psphot/src/psphot.c =================================================================== --- trunk/psphot/src/psphot.c (revision 20356) +++ trunk/psphot/src/psphot.c (revision 20411) @@ -10,6 +10,5 @@ psTimerStart ("complete"); pmErrorRegister(); // register psModule's error codes/messages - psphotErrorRegister(); // register our error codes/messages - psphotModelClassInit (); // load implementation-specific models + psphotInit(); // load command-line arguments, options, and system config data Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 20356) +++ trunk/psphot/src/psphot.h (revision 20411) @@ -13,4 +13,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 @@ -20,4 +28,5 @@ bool psphotModelTest (pmConfig *config, const pmFPAview *view, psMetadata *recipe); +bool psphotInit (); bool psphotReadout (pmConfig *config, const pmFPAview *view); bool psphotReadoutFindPSF(pmConfig *config, const pmFPAview *view); @@ -41,5 +50,6 @@ bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources); bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final); -bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf); +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); @@ -51,4 +61,5 @@ bool psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe); bool psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe); +bool psphotSetThreads (); // used by psphotFindDetections Index: trunk/psphot/src/psphotArguments.c =================================================================== --- trunk/psphot/src/psphotArguments.c (revision 20356) +++ trunk/psphot/src/psphotArguments.c (revision 20411) @@ -29,4 +29,17 @@ // these options override the PSPHOT recipe values loaded from recipe files psMetadata *options = pmConfigRecipeOptions (config, PSPHOT_RECIPE); + + // Number of threads + if ((N = psArgumentGet(argc, argv, "-threads"))) { + psArgumentRemove(N, &argc, argv); + int nThreads = atoi(argv[N]); + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "NTHREADS", 0, "number of psphot threads", nThreads); + psArgumentRemove(N, &argc, argv); + + // 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); + } + psphotSetThreads(); // run the test model (requires X,Y coordinate) Index: trunk/psphot/src/psphotCleanup.c =================================================================== --- trunk/psphot/src/psphotCleanup.c (revision 20356) +++ trunk/psphot/src/psphotCleanup.c (revision 20411) @@ -5,4 +5,5 @@ psFree (config); + psThreadPoolFinalize (); psTimerStop (); psMemCheckCorruption (stderr, true); Index: trunk/psphot/src/psphotGuessModels.c =================================================================== --- trunk/psphot/src/psphotGuessModels.c (revision 20356) +++ trunk/psphot/src/psphotGuessModels.c (revision 20411) @@ -1,3 +1,7 @@ # include "psphotInternal.h" + +bool psphotChooseCellSizes (int *Cx, int *Cy, pmReadout *readout, int nThreads); +psphotGuessModelForRegionArgs *psphotGuessModelForRegionArgsAlloc(); +bool psphotGuessModelForRegion (psphotGuessModelForRegionArgs *args); // A guess for when the moments aren't available @@ -16,10 +20,22 @@ } +# define NFILL 4 + // construct an initial PSF model for each object -bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { +bool psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) { bool status; psTimerStart ("psphot"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + assert (recipe); + + // int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + // if (!status) { + // nThreads = 0; + // } + int nThreads = 0; // bit-masks to test for good/bad pixels @@ -37,4 +53,133 @@ 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: + 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)) { + 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); + } + } + } + + psLogMsg ("psphot.models", 4, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot")); + 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; + + int nSrc = 0; + for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; @@ -43,5 +188,13 @@ if (source->type == PM_SOURCE_TYPE_DEFECT) continue; 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 ++; + // XXX if a source is faint, it will not have moments measured. // it must be modelled as a PSF. In this case, we need to use @@ -53,5 +206,5 @@ } else { // use the source moments, etc to guess basic model parameters - modelEXT = pmSourceModelGuess (source, psf->type); + modelEXT = pmSourceModelGuess (source, psf->type); // ALLOC if (!modelEXT) { modelEXT = wildGuess(source, psf); @@ -68,5 +221,5 @@ // set PSF parameters for this model (apply 2D shape model) - pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); + pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); // ALLOC if (modelPSF == NULL) { psError(PSPHOT_ERR_PSF, false, @@ -92,4 +245,5 @@ // XXX need to define the guess flux? // set the fit radius based on the object flux limit and the model + // this function affects the mask pixels psphotCheckRadiusPSF (readout, source, modelPSF, markVal); @@ -99,8 +253,35 @@ pmSourceCacheModel (source, maskVal); + } - psLogMsg ("psphot.models", 4, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot")); + + // fprintf (stderr, "%d for region %lf,%lf - %lf,%lf\n", nSrc, region->x0, region->y0, region->x1, region->y1); return true; } -// XXX do we always know which model is supposed to be used? +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; +} Index: trunk/psphot/src/psphotReadout.c =================================================================== --- trunk/psphot/src/psphotReadout.c (revision 20356) +++ trunk/psphot/src/psphotReadout.c (revision 20411) @@ -1,3 +1,12 @@ # include "psphotInternal.h" + +// this should be called by every program that links against libpsphot +bool psphotInit () { + + psphotErrorRegister(); // register our error codes/messages + psphotModelClassInit (); // load implementation-specific models + psphotSetThreads (); + return true; +} bool psphotReadout(pmConfig *config, const pmFPAview *view) { @@ -149,5 +158,5 @@ // construct an initial model for each object - psphotGuessModels (readout, sources, recipe, psf); + psphotGuessModels (config, readout, sources, psf); // XXX test output of models @@ -207,5 +216,5 @@ // create full input models - psphotGuessModels (readout, newSources, recipe, psf); + psphotGuessModels (config, readout, newSources, psf); // replace all sources so fit below applies to all at once Index: trunk/psphot/src/psphotSetThreads.c =================================================================== --- trunk/psphot/src/psphotSetThreads.c (revision 20411) +++ trunk/psphot/src/psphotSetThreads.c (revision 20411) @@ -0,0 +1,22 @@ +# 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 () { + + psThreadTask *task = NULL; + + task = psThreadTaskAlloc("PSPHOT_GUESS_MODEL", 1); + task->function = &psphotThread_psphotGuessModel; + psThreadTaskAdd(task); + psFree(task); + + return true; +}