Index: /trunk/ppMerge/src/Makefile.am =================================================================== --- /trunk/ppMerge/src/Makefile.am (revision 18838) +++ /trunk/ppMerge/src/Makefile.am (revision 18839) @@ -10,11 +10,11 @@ ppMergeFiles.c \ ppMergeScaleZero.c \ - ppMergeLoop.c \ ppMergeFileGroup.c \ + ppMergeReadChunk.c \ + ppMergeLoop_Threaded.c \ + ppMergeSetThreads.c \ ppMergeMask.c -# ppMergeLoop_Threaded.c -# ppMergeThreadLauncher.c - +# ppMergeLoop.c noinst_HEADERS = \ Index: /trunk/ppMerge/src/ppMerge.c =================================================================== --- /trunk/ppMerge/src/ppMerge.c (revision 18838) +++ /trunk/ppMerge/src/ppMerge.c (revision 18839) @@ -10,5 +10,4 @@ { psLibInit(NULL); - psMemSetThreadSafety(false); psTimerStart(TIMERNAME); Index: /trunk/ppMerge/src/ppMerge.h =================================================================== --- /trunk/ppMerge/src/ppMerge.h (revision 18838) +++ /trunk/ppMerge/src/ppMerge.h (revision 18839) @@ -15,5 +15,5 @@ #define TIMERNAME "ppMerge" // Name for timer #define PPMERGE_RECIPE "PPMERGE" // Recipe name -#define THREADED 0 +#define THREADED 1 // Type of frame to merge @@ -39,4 +39,6 @@ bool read; bool busy; + int firstScan; + int lastScan; } ppMergeFileGroup; @@ -100,6 +102,8 @@ ppMergeFileGroup *ppMergeFileGroupAlloc(); -ppMergeFileGroup *ppMergeReadChunk (psArray *fileGroups, pmConfig *config, int numChunk); +ppMergeFileGroup *ppMergeReadChunk (bool *status, psArray *fileGroups, pmConfig *config, int numChunk); void *ppMergeThreadLauncher (void *data); +bool ppMergeSetThreads (); + #endif Index: /trunk/ppMerge/src/ppMergeArguments.c =================================================================== --- /trunk/ppMerge/src/ppMergeArguments.c (revision 18838) +++ /trunk/ppMerge/src/ppMergeArguments.c (revision 18839) @@ -171,5 +171,4 @@ } -# if (THREADED) // Number of threads if ((argnum = psArgumentGet(argc, argv, "-threads"))) { @@ -181,7 +180,7 @@ // 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, &ppMergeThreadLauncher); - } -# endif + psThreadPoolInit (nThreads); + } + ppMergeSetThreads(); if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 3) { Index: /trunk/ppMerge/src/ppMergeLoop_Threaded.c =================================================================== --- /trunk/ppMerge/src/ppMergeLoop_Threaded.c (revision 18838) +++ /trunk/ppMerge/src/ppMergeLoop_Threaded.c (revision 18839) @@ -36,5 +36,4 @@ // General combination parameters - int rows = psMetadataLookupS32(NULL, arguments, "ROWS"); // Number of rows to read per chunk int iter = psMetadataLookupS32(NULL, arguments, "ITER"); // Number of rejection iterations float rej = psMetadataLookupF32(NULL, arguments, "REJ"); // Rejection level @@ -175,57 +174,72 @@ fileGroup->read = false; fileGroup->busy = false; + fileGroup->lastScan = 0; + fileGroup->firstScan = 0; fileGroups->data[i] = fileGroup; } + // call the init functions + switch (type) { + case PPMERGE_TYPE_BIAS: + case PPMERGE_TYPE_FLAT: + case PPMERGE_TYPE_FRINGE: + psAssert (fileGroups->n > 0, "no valid file groups defined"); + ppMergeFileGroup *fileGroup = fileGroups->data[0]; + if (!pmReadoutCombinePrepare(outRO, fileGroup->readouts, combination)) { + goto ERROR; + } + break; + + default: + fprintf (stderr, "not yet ready"); + goto ERROR; + } + // Read input data by chunks + // psTimerStart ("ppMergeLoop"); for (int numChunk = 0; true; numChunk++) { - + + bool status = false; ppMergeFileGroup *fileGroup = ppMergeReadChunk (&status, fileGroups, config, numChunk); if (!status) goto ERROR; if (!fileGroup) break; + psThreadJob *job = NULL; + switch (type) { case PPMERGE_TYPE_SHUTTER: - if (nThreads) { - // allocate a job - psThreadJob *job = psThreadJobAlloc ("PPMERGE_SHUTTER_CORRECTION", 0); - - // construct the arguments for this job - psArrayAdd (job->args, 1, outRO); - psArrayAdd (job->args, 1, fileGroup); - psArrayAdd (job->args, 1, psScalarAlloc(shutterRef, PS_TYPE_F32)) - psArrayAdd (job->args, 1, shutters->data[cellNum]); - psArrayAdd (job->args, 1, psScalarAlloc(iter, PS_TYPE_S32)); - psArrayAdd (job->args, 1, psScalarAlloc(rej, PS_TYPE_F32)); - psArrayAdd (job->args, 1, psScalarAlloc(maskVal, PS_TYPE_U8)); - - psThreadJobAddPending (job); - } else { - if (!pmShutterCorrectionGenerate(outRO, NULL, fileGroup->readouts, shutterRef, shutters->data[cellNum], iter, rej, maskVal)) { - goto ERROR; - } - fileGroup->busy = false; + // allocate a job + job = psThreadJobAlloc ("PPMERGE_SHUTTER_CORRECTION"); + + // construct the arguments for this job + psArrayAdd (job->args, 1, outRO); + psArrayAdd (job->args, 1, fileGroup); + psArrayAdd (job->args, 1, psScalarAlloc(shutterRef, PS_TYPE_F32)); + psArrayAdd (job->args, 1, shutters->data[cellNum]); + psArrayAdd (job->args, 1, psScalarAlloc(iter, PS_TYPE_S32)); + psArrayAdd (job->args, 1, psScalarAlloc(rej, PS_TYPE_F32)); + psArrayAdd (job->args, 1, psScalarAlloc(maskVal, PS_TYPE_U8)); + + // call: pmShutterCorrectionGenerate(outRO, NULL, fileGroup->readouts, shutterRef, shutters->data[cellNum], iter, rej, maskVal) + if (!psThreadJobAddPending (job)) { + goto ERROR; } break; case PPMERGE_TYPE_DARK: - if (nThreads) { - // allocate a job - psThreadJob *job = psThreadJobAlloc ("PPMERGE_DARK_COMBINE", 0); - - // construct the arguments for this job - psArrayAdd (job->args, 1, outCell); - psArrayAdd (job->args, 1, fileGroup); - psArrayAdd (job->args, 1, darkOrdinates); - psArrayAdd (job->args, 1, darkNorm); - psArrayAdd (job->args, 1, psScalarAlloc(iter, PS_TYPE_S32)); - psArrayAdd (job->args, 1, psScalarAlloc(rej, PS_TYPE_F32)); - psArrayAdd (job->args, 1, psScalarAlloc(maskVal, PS_TYPE_U8)); - - psThreadJobAddPending (job); - } else { - if (!pmDarkCombine(outCell, fileGroup->readouts, darkOrdinates, darkNorm, iter, rej, maskVal)) { - goto ERROR; - } - fileGroup->busy = false; + // allocate a job + job = psThreadJobAlloc ("PPMERGE_DARK_COMBINE"); + + // construct the arguments for this job + psArrayAdd (job->args, 1, outCell); + psArrayAdd (job->args, 1, fileGroup); + psArrayAdd (job->args, 1, darkOrdinates); + psArrayAdd (job->args, 1, darkNorm); + psArrayAdd (job->args, 1, psScalarAlloc(iter, PS_TYPE_S32)); + psArrayAdd (job->args, 1, psScalarAlloc(rej, PS_TYPE_F32)); + psArrayAdd (job->args, 1, psScalarAlloc(maskVal, PS_TYPE_U8)); + + // call: pmDarkCombine(outCell, fileGroup->readouts, darkOrdinates, darkNorm, iter, rej, maskVal); + if (!psThreadJobAddPending (job)) { + goto ERROR; } break; @@ -233,21 +247,17 @@ case PPMERGE_TYPE_FLAT: case PPMERGE_TYPE_FRINGE: - if (nThreads) { - // allocate a job - psThreadJob *job = psThreadJobAlloc ("PPMERGE_READOUT_COMBINE", 0); - - // construct the arguments for this job - psArrayAdd (job->args, 1, outRO); - psArrayAdd (job->args, 1, fileGroup); - psArrayAdd (job->args, 1, zeros); - psArrayAdd (job->args, 1, scales); - psArrayAdd (job->args, 1, combination); - - psThreadJobAddPending (job); - } else { - if (!pmReadoutCombine(outRO, fileGroup->readouts, zeros, scales, combination)) { - goto ERROR; - } - fileGroup->busy = false; + // allocate a job + job = psThreadJobAlloc ("PPMERGE_READOUT_COMBINE"); + + // construct the arguments for this job + psArrayAdd (job->args, 1, outRO); + psArrayAdd (job->args, 1, fileGroup); + psArrayAdd (job->args, 1, zeros); + psArrayAdd (job->args, 1, scales); + psArrayAdd (job->args, 1, combination); + + // call: pmReadoutCombine(outRO, fileGroup->readouts, zeros, scales, combination); + if (!psThreadJobAddPending (job)) { + goto ERROR; } break; @@ -258,24 +268,23 @@ // wait for the threads to finish and manage results - if (nThreads) { - // wait here for the threaded jobs to finish - if (!psThreadPoolWait ()) { - psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); - return false; - } - fprintf (stderr, "success for threaded jobs\n"); - - // we don't care about the results, just dump the done queue jobs - psThreadJob *job = NULL; - while ((job = psThreadJobGetDone()) != NULL) { - psFree (job); - } + if (!psThreadPoolWait ()) { + psError(PS_ERR_UNKNOWN, false, "Unable to combine images."); + return false; } + + // we don't care about the results, just dump the done queue jobs + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + psFree (job); + } + + psFree(fileGroups); // Get list of cells for concepts averaging psList *inCells = psListAlloc(NULL); // List of cells for (int i = 0; i < numFiles; i++) { - pmReadout *readout = readouts->data[i]; // Readout of interest - psListAdd(inCells, PS_LIST_TAIL, readout->parent); + pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); + pmCell *inCell = pmFPAviewThisCell(view, input->fpa); // Input cell + psListAdd(inCells, PS_LIST_TAIL, inCell); } if (!pmConceptsAverageCells(outCell, inCells, NULL, NULL, true)) { @@ -286,6 +295,5 @@ } psFree(inCells); - - psFree(fileGroups); + // fprintf (stdout, "done ppMergeLoop for cell : %f\n", psTimerMark ("ppMergeLoop")); // Plug supplementary images into their own FPAs Index: /trunk/ppMerge/src/ppMergeReadChunk.c =================================================================== --- /trunk/ppMerge/src/ppMergeReadChunk.c (revision 18838) +++ /trunk/ppMerge/src/ppMergeReadChunk.c (revision 18839) @@ -1,11 +1,14 @@ # include "ppMerge.h" -ppMergeFileGroup *ppMergeReadChunk (psArray *fileGroups, pmConfig *config, int numChunk) { +ppMergeFileGroup *ppMergeReadChunk (bool *status, psArray *fileGroups, pmConfig *config, int numChunk) { + + *status = true; + + bool mdok; + bool haveMasks = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS"); // Do we have masks? + bool haveWeights = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS"); // Do we have weights? + int rows = psMetadataLookupS32(NULL, config->arguments, "ROWS"); // Number of rows to read per chunk // select an available fileGroup - - bool haveMasks = psMetadataLookupBool(&mdok, arguments, "INPUTS.MASKS"); // Do we have masks? - bool haveWeights = psMetadataLookupBool(&mdok, arguments, "INPUTS.WEIGHTS"); // Do we have weights? - while (1) { // check for any fileGroups which can read data @@ -14,5 +17,16 @@ if (fileGroup->read) continue; + // find max last scan so far + int lastScan = 0; + for (int i = 0; i < fileGroups->n; i++) { + ppMergeFileGroup *fileGroup = fileGroups->data[i]; + lastScan = PS_MAX (fileGroup->lastScan, lastScan); + } + fileGroup->firstScan = lastScan; + fileGroup->lastScan = lastScan + rows; + psArray *readouts = fileGroup->readouts; + + psTimerStart ("ppMergeReadChunk"); psTrace("ppStack", 2, "Reading data for chunk %d into fileGroup %d....n", numChunk, j); @@ -20,27 +34,46 @@ pmReadout *inRO = readouts->data[i]; // Input readout + // override the recorded last scan + inRO->thisImageScan = fileGroup->firstScan; + inRO->thisWeightScan = fileGroup->firstScan; + inRO->thisMaskScan = fileGroup->firstScan; + // Read a chunk from a file pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); - if (!pmReadoutReadChunk(inRO, file->fits, 0, rows, 0, config)) { - psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPMERGE.INPUT %d", numChunk, i); - return NULL; - } - if (haveMasks) { - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.MASK", i); - if (!pmReadoutReadChunkMask(inRO, file->fits, 0, rows, 0, config)) { - psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPMERGE.INPUT.MASK %d", numChunk, NAME, i); + bool keepReading = false; + if (pmReadoutMore(inRO, file->fits, 0, rows, config)) { + keepReading = true; + if (!pmReadoutReadChunk(inRO, file->fits, 0, rows, 0, config)) { + psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPMERGE.INPUT %d", numChunk, i); + *status = false; return NULL; } } - if (haveWeights) { - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.WEIGHT", i); - if (!pmReadoutReadChunkWeight(inRO, file->fits, 0, rows, 0, config)) { - psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPMERGE.INPUT.WEIGHT %d", numChunk, NAME, i); + if (haveMasks && pmReadoutMoreMask(inRO, file->fits, 0, rows, config)) { + keepReading = true; + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.MASK", i); + if (!pmReadoutReadChunkMask(inRO, file->fits, 0, rows, 0, config)) { + psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPMERGE.INPUT.MASK %d", numChunk, i); + *status = false; return NULL; } } + + if (haveWeights && pmReadoutMoreWeight(inRO, file->fits, 0, rows, config)) { + keepReading = true; + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.WEIGHT", i); + if (!pmReadoutReadChunkWeight(inRO, file->fits, 0, rows, 0, config)) { + psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPMERGE.INPUT.WEIGHT %d", numChunk, i); + *status = false; + return NULL; + } + } + if (!keepReading) { + return NULL; + } } + fileGroup->read = fileGroup->busy = true; return fileGroup; @@ -48,32 +81,11 @@ // check for any fileGroups which are done processing - bool wait = true; - bool more = true; + bool wait = false; for (int j = 0; j < fileGroups->n; j++) { ppMergeFileGroup *fileGroup = fileGroups->data[j]; if (!fileGroup->read || fileGroup->busy) continue; - - wait = false; - psArray *readouts = fileGroup->readouts; - // any more data to be read? - for (int i = 0; i < readouts->n && more; i++) { - pmReadout *inRO = readouts->data[i]; - - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); - more &= pmReadoutMore(inRO, file->fits, 0, rows, config); - - if (haveMasks) { - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.MASK", i); - more &= pmReadoutMoreMask(inRO, file->fits, 0, rows, config); - } - if (haveWeights) { - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.WEIGHT", i); - more &= pmReadoutMoreWeight(inRO, file->fits, 0, rows, config); - } - } fileGroup->read = false; + wait = true; } - if (!more) return NULL; - if (wait) usleep (10000); } Index: /trunk/ppMerge/src/ppMergeSetThreads.c =================================================================== --- /trunk/ppMerge/src/ppMergeSetThreads.c (revision 18839) +++ /trunk/ppMerge/src/ppMergeSetThreads.c (revision 18839) @@ -0,0 +1,70 @@ +# include "ppMerge.h" + +// "PPMERGE_READOUT_COMBINE", 5 +bool ppMergeThread_pmReadoutCombine (psThreadJob *job) { + + pmReadout *output = job->args->data[0]; + ppMergeFileGroup *fileGroup = job->args->data[1]; + psVector *zero = job->args->data[2]; + psVector *scale = job->args->data[3]; + pmCombineParams *params = job->args->data[4]; + + bool status = pmReadoutCombine (output, fileGroup->readouts, zero, scale, params); + + // after we are done, tell the I/O system that this file group is done + fileGroup->busy = false; + return status; +} + +bool ppMergeThread_pmDarkCombine (psThreadJob *job) { + + pmCell *outCell = job->args->data[0]; + ppMergeFileGroup *fileGroup = job->args->data[1]; + psArray *darkOrdinates = job->args->data[2]; + psString darkNorm = job->args->data[3]; + psScalar *iter = job->args->data[4]; + psScalar *rej = job->args->data[5]; + psScalar *maskVal = job->args->data[6]; + + bool status = pmDarkCombine(outCell, fileGroup->readouts, darkOrdinates, darkNorm, iter->data.S32, rej->data.F32, maskVal->data.U8); + + // after we are done, tell the I/O system that this file group is done + fileGroup->busy = false; + return status; +} + +bool ppMergeThread_pmShuttderCorrectionGenerate (psThreadJob *job) { + + pmReadout *output = job->args->data[0]; + ppMergeFileGroup *fileGroup = job->args->data[1]; + psScalar *shutterRef = job->args->data[2]; + pmShutterCorrectionData *data = job->args->data[3]; + psScalar *iter = job->args->data[4]; + psScalar *rej = job->args->data[5]; + psScalar *maskVal = job->args->data[6]; + + bool status = pmShutterCorrectionGenerate(output, NULL, fileGroup->readouts, shutterRef->data.F32, data, iter->data.S32, rej->data.F32, maskVal->data.U8); + + // after we are done, tell the I/O system that this file group is done + fileGroup->busy = false; + return status; +} + +bool ppMergeSetThreads () { + + psThreadTask *task = NULL; + + task = psThreadTaskAlloc ("PPMERGE_READOUT_COMBINE", 5); + task->function = &ppMergeThread_pmReadoutCombine; + psThreadTaskAdd (task); + + task = psThreadTaskAlloc ("PPMERGE_DARK_COMBINE", 7); + task->function = &ppMergeThread_pmReadoutCombine; + psThreadTaskAdd (task); + + task = psThreadTaskAlloc ("PPMERGE_SHUTTER_CORRECTION", 7); + task->function = &ppMergeThread_pmReadoutCombine; + psThreadTaskAdd (task); + + return true; +} Index: /trunk/ppMerge/src/ppMergeThreadLauncher.c =================================================================== --- /trunk/ppMerge/src/ppMergeThreadLauncher.c (revision 18838) +++ /trunk/ppMerge/src/ppMergeThreadLauncher.c (revision 18839) @@ -41,4 +41,5 @@ self->fault = true; } + // we do not have to lock here because this transition is not tied to the job queue fileGroup->busy = false; Index: /trunk/pswarp/src/Makefile.am =================================================================== --- /trunk/pswarp/src/Makefile.am (revision 18838) +++ /trunk/pswarp/src/Makefile.am (revision 18839) @@ -16,5 +16,8 @@ pswarpPixelFraction.c \ pswarpSetMaskBits.c \ - pswarpTransformReadout_Opt.c \ + pswarpSetThreads.c \ + pswarpTransformReadout.c \ + pswarpTransformSources.c \ + pswarpTransformTile.c \ pswarpVersion.c Index: /trunk/pswarp/src/pswarp.h =================================================================== --- /trunk/pswarp/src/pswarp.h (revision 18838) +++ /trunk/pswarp/src/pswarp.h (revision 18839) @@ -10,4 +10,6 @@ #include #include + +#define THREADED 1 #include "pswarpErrorCodes.h" @@ -33,6 +35,25 @@ int nXpts, nYpts; // number of x,y samples in the grid int nXpix, nYpix; // x,y spacing in src image pixels of grid samples - int xMin, yMin; // coordinate of first grid sample + double xMin, yMin; // coordinate of first grid sample } pswarpMapGrid; + +typedef struct { + // values which are common to all tiles + pmReadout *input; + pmReadout *output; + pswarpMapGrid *grid; + psImageInterpolateOptions *interp; + psImage *region; + + // input values for this tile + int gridX; + int gridY; + + // output values for this tile + long goodPixels; +} pswarpTransformTileArgs; + +pswarpTransformTileArgs *pswarpTransformTileArgsAlloc(); +bool pswarpTransformTile (pswarpTransformTileArgs *args); pmConfig *pswarpArguments (int argc, char **argv); @@ -44,5 +65,5 @@ void pswarpCleanup (pmConfig *config); bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config); -bool pswarpTransformReadout_Opt (pmReadout *output, pmReadout *input, pmConfig *config); +bool pswarpTransformSources(pmReadout *output, pmReadout *input, pmConfig *config); bool pswarpMatchRange (int *minX, int *minY, int *maxX, int *maxY, pmReadout *dest, pmReadout *src); @@ -53,4 +74,5 @@ pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix); bool pswarpMapGridSetGrid (pswarpMapGrid *grid, int ix, int iy, int *gridX, int *gridY); +bool pswarpMapGridCoordRange (pswarpMapGrid *grid, int gridX, int gridY, psPlane *min, psPlane *max); int pswarpMapGridNextGrid_X (pswarpMapGrid *grid, int gridX); int pswarpMapGridNextGrid_Y (pswarpMapGrid *grid, int gridY); @@ -68,2 +90,4 @@ ); +// define threads for this program +bool pswarpSetThreads (); Index: /trunk/pswarp/src/pswarpArguments.c =================================================================== --- /trunk/pswarp/src/pswarpArguments.c (revision 18838) +++ /trunk/pswarp/src/pswarpArguments.c (revision 18839) @@ -40,4 +40,17 @@ psArgumentRemove(N, &argc, argv); } + + // 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 warp 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); + } + pswarpSetThreads (); // PSF determination? Index: /trunk/pswarp/src/pswarpCleanup.c =================================================================== --- /trunk/pswarp/src/pswarpCleanup.c (revision 18838) +++ /trunk/pswarp/src/pswarpCleanup.c (revision 18839) @@ -9,4 +9,5 @@ pmModelClassCleanup (); psTimeFinalize (); + psThreadPoolFinalize (); pmConceptsDone (); pmConfigDone (); Index: /trunk/pswarp/src/pswarpLoop.c =================================================================== --- /trunk/pswarp/src/pswarpLoop.c (revision 18838) +++ /trunk/pswarp/src/pswarpLoop.c (revision 18839) @@ -241,5 +241,5 @@ } - pswarpTransformReadout_Opt(output, readout, config); + pswarpTransformReadout(output, readout, config); pmFPAfileIOChecks(config, view, PM_FPA_AFTER); Index: /trunk/pswarp/src/pswarpMapGrid.c =================================================================== --- /trunk/pswarp/src/pswarpMapGrid.c (revision 18838) +++ /trunk/pswarp/src/pswarpMapGrid.c (revision 18839) @@ -4,5 +4,5 @@ // coordinates (src) to destination coordinates (dest). we construct a grid with superpixel // spacing of nXpix, nYpix. The transformation for each grid cell is valid for the superpixel. -// The grid over-fills the source image so ever source image pixel is guaranteed to have a map. +// The grid over-fills the source image so every source image pixel is guaranteed to have a map. pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix) { @@ -11,6 +11,6 @@ // start counting from the center of the superpixels - int xMin = 0.5*nXpix; - int yMin = 0.5*nYpix; + double xMin = 0.5*nXpix; + double yMin = 0.5*nYpix; // the map is defined for coordinates in the image parent frame. @@ -42,8 +42,22 @@ // set the grid coordinate (gridX,gridY) for the given source image coordinate (ix,iy) +// XXX return true if the result is on the src image, false otherwise (???) bool pswarpMapGridSetGrid (pswarpMapGrid *grid, int ix, int iy, int *gridX, int *gridY) { *gridX = 0.5 + (ix - grid->xMin) / (double) grid->nXpix; *gridY = 0.5 + (iy - grid->yMin) / (double) grid->nYpix; + + return true; +} + +// given the specified grid coordinate (gridX, gridY), return the min and max coordinates for the tile +bool pswarpMapGridCoordRange (pswarpMapGrid *grid, int gridX, int gridY, psPlane *min, psPlane *max) { + + min->x = (gridX - 0.5)*grid->nXpix + grid->xMin; + min->y = (gridY - 0.5)*grid->nYpix + grid->yMin; + + max->x = min->x + grid->nXpix; + max->y = min->y + grid->nYpix; + return true; } Index: /trunk/pswarp/src/pswarpSetThreads.c =================================================================== --- /trunk/pswarp/src/pswarpSetThreads.c (revision 18839) +++ /trunk/pswarp/src/pswarpSetThreads.c (revision 18839) @@ -0,0 +1,21 @@ +# include "pswarp.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 pswarpThread_pswarpTransformTile (psThreadJob *job) { + + pswarpTransformTileArgs *args = job->args->data[0]; + bool status = pswarpTransformTile (args); + return status; +} + +bool pswarpSetThreads () { + + psThreadTask *task = NULL; + + task = psThreadTaskAlloc ("PSWARP_TRANSFORM_TILE", 1); + task->function = &pswarpThread_pswarpTransformTile; + psThreadTaskAdd (task); + + return true; +} Index: /trunk/pswarp/src/pswarpThreadLauncher.c =================================================================== --- /trunk/pswarp/src/pswarpThreadLauncher.c (revision 18839) +++ /trunk/pswarp/src/pswarpThreadLauncher.c (revision 18839) @@ -0,0 +1,39 @@ +# include "pswarp.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 +void *pswarpThreadLauncher (void *data) { + + psThread *self = data; + psThreadJob *job = NULL; + + while (1) { + + // if we get an error, just wait until we are cleared or killed + while (self->fault) { + usleep (10000); + } + + // request a new job, if there are none available, sleep a bit + // we have to lock here so the job queue cannot be empty yet no threads busy + psThreadLock(); + while ((job = psThreadJobGetPending ()) == NULL) { + psThreadUnlock(); + usleep (10000); + } + self->busy = true; + psThreadUnlock(); + + // list all allowed job types here + if (!strcmp (job->type, "PSWARP_TRANSFORM_TILE")) { + pswarpTransformTileArgs *args = job->args->data[0]; + bool status = pswarpTransformTile (args); + if (!status) { + self->fault = true; + } + // we do not have to lock here because this transition is not tied to the job queue + self->busy = false; + continue; + } + } +} Index: /trunk/pswarp/src/pswarpThreads.c =================================================================== --- /trunk/pswarp/src/pswarpThreads.c (revision 18839) +++ /trunk/pswarp/src/pswarpThreads.c (revision 18839) @@ -0,0 +1,52 @@ +# include "pswarp.h" + +typedef enum { + PSWARP_TRANSFORM_TILE +} pswarpJobType; + +typedef struct { + pswarpJobType type; + void *opts; +} pswarpJob; + +static psArray *jobs = NULL; +static pthread_t *threads = NULL; + +bool pswarpThreadsAddJob (pwarpJob *job) { + + if (jobs == NULL) { + jobs = psArrayAlloc (16); + } + + psArrayAdd (jobs, job); + return true; +} + +bool pswarpThreadsLaunchJobs () { + + while () { + while ((job = psListGetAndRemove (jobs, PS_LIST_HEAD)) == NULL) { + usleep (50000); + } + + switch (job->type) { + case PSWARP_TRANSFORM_TILE: + status = pswarpTransformTile ((pswarpTransformTileOpts *)job->opts); + // send a message somewhere if the job fails + break; + } + } +} + +// create a pool of Nthreads +bool pswarpThreadsCreate (int nThreads) { + + threads = (pthread_t *) psAlloc (sizeof(pthread_t)); + + for (int i = 0; i < nThreads; i++) { + pthread_create (&threads[i], NULL, &function, NULL); + } + + return true; +} + Index: /trunk/pswarp/src/pswarpTransformReadout.c =================================================================== --- /trunk/pswarp/src/pswarpTransformReadout.c (revision 18838) +++ /trunk/pswarp/src/pswarpTransformReadout.c (revision 18839) @@ -1,78 +1,150 @@ # include "pswarp.h" -bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config) { - +// NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT +bool pswarpTransformReadout(pmReadout *output, pmReadout *input, pmConfig *config) +{ // XXX this implementation currently ignores the use of the region psImage *region = NULL; - pmCell *cell = NULL; - // select the current recipe - // psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE); + psTimerStart ("warp"); - int outNx = output->image->numCols; - int outNy = output->image->numRows; + // Get warp parameters + bool mdok; + int nGridX = psMetadataLookupS32(NULL, config->arguments, "GRID.NX"); + int nGridY = psMetadataLookupS32(NULL, config->arguments, "GRID.NY"); + psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.MODE"); - psPlane *inPix = psPlaneAlloc(); // Coordinates on the input detector - psPlane *outPix = psPlaneAlloc(); // Coordinates on the output detector + // load the recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE); + psAssert (recipe, "missing recipe %s", PSWARP_RECIPE); - psPlane *FP = psPlaneAlloc(); // Coordinates on the focal plane - psPlane *TP = psPlaneAlloc(); // Coordinates on the tangent plane - psSphere *sky = psSphereAlloc(); // Coordinates on the sky + // output mask bits + psMaskType maskIn = psMetadataLookupU8(&mdok, recipe, "MASK.INPUT"); + psMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config); + psMaskType maskBad = pmConfigMaskGet("BAD.WARP", config); + psAssert (mdok, "MASK.INPUT was not defined"); - cell = input->parent; - pmChip *chipInput = cell->parent; - pmFPA *fpaInput = chipInput->parent; + int nThreads = psMetadataLookupS32 (&mdok, config->arguments, "NTHREADS"); + if (!mdok) nThreads = 0; - cell = output->parent; - pmChip *chipOutput = cell->parent; - pmFPA *fpaOutput = chipOutput->parent; + // Flux fraction for "poor" + float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC"); - psF32 **outData = output->image->data.F32; - psImage *inImage = input->image; + // pswarpMapGridFromImage builds a set of locally-linear maps which convert the + // output coordinates to input coordinates + pswarpMapGrid *grid = pswarpMapGridFromImage (input, output, nGridX, nGridY); - psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR, inImage, - NULL, NULL, 0, NAN, NAN, 0, 0, 0.0); + // XXX optionally modify the grid based on this result and force the maxError < XXX + double maxError = pswarpMapGridMaxError (grid); + psLogMsg ("pswarp", 3, "maximum error using this grid sampling: %f\n", maxError); - // Iterate over the output image pixels - for (int y = 0; y < outNy; y++) { - for (int x = 0; x < outNx; x++) { - // Only transform those pixels requested - if (region && region->data.U8[y][x]) continue; + // Interpolation options : move these from the arguments to explicit assignments + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(interpolationMode, input->image, input->weight, input->mask, maskIn, NAN, NAN, maskBad, maskPoor, poorFrac); - // XXX double check this 1/2 pixel offset - outPix->x = (double)x + 0.5; - outPix->y = (double)y + 0.5; + if (input->weight && !output->weight) { + output->weight = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_F32); + psImageInit(output->weight, NAN); + } + if ((input->mask || maskPoor || maskBad) && !output->mask) { + output->mask = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_MASK); + psImageInit(output->mask, maskBad); + } - psPlaneTransformApply(FP, chipOutput->toFPA, outPix); - psPlaneTransformApply (TP, fpaOutput->toTPA, FP); - psDeproject (sky, TP, fpaOutput->toSky); + // total number of good pixels across all tiles (summed below) + int goodPixels = 0; - psProject (TP, sky, fpaInput->toSky); - psPlaneTransformApply (FP, fpaInput->fromTPA, TP); - psPlaneTransformApply (inPix, chipInput->fromFPA, FP); + // create jobs and supply them to the threads + for (int gridY = 0; gridY < grid->nYpts; gridY++) { + for (int gridX = 0; gridX < grid->nXpts; gridX++) { - // XXX get interpolation method from the recipe - double value; - if (!psImageInterpolate(&value, NULL, NULL, inPix->x, inPix->y, interp)) { - psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); - psFree(interp); - psFree(inPix); - psFree(outPix); - psFree(FP); - psFree(TP); - psFree(sky); - return false; - } + pswarpTransformTileArgs *args = pswarpTransformTileArgsAlloc(); - outData[y][x] = value; - // modify zero and scale? + // these items are just views to the data; they are not freed with args + args->input = input; + args->output = output; + args->grid = grid; + args->interp = interp; + args->region = region; + + args->gridX = gridX; + args->gridY = gridY; + args->goodPixels = 0; + + // allocate a job + psThreadJob *job = psThreadJobAlloc ("PSWARP_TRANSFORM_TILE"); + + // construct the arguments for this job + // job is pswarpTransformTile (gridX, gridY); + psArrayAdd (job->args, 1, args); + // fprintf (stderr, "adding job %d,%d, Nargs: %ld\n", gridX, gridY, job->args->n); + + // call: pswarpTransformTile (args); + if (!psThreadJobAddPending (job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to warp image."); + return false; + } + psFree (args); + } + } + + // wait for the threads to finish and manage results + // wait here for the threaded jobs to finish + if (!psThreadPoolWait ()) { + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); + return false; + } + fprintf (stderr, "success for threaded jobs\n"); + + // each job records its own goodPixel values; sum them here + // 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 { + pswarpTransformTileArgs *args = job->args->data[0]; + // fprintf (stderr, "finished job %d,%d, Nargs: %ld\n", args->gridX, args->gridY, job->args->n); + goodPixels += args->goodPixels; + } + psFree (job); + } + psFree(interp); + psFree(grid); + + + // Store the variance factor and number of good pixels + if (goodPixels > 0) { + // Variance factor: large factor --> small scale + float varFactor = psImageInterpolateVarianceFactor(input->image->numCols / 2.0 + input->image->col0, input->image->numRows / 2.0 + input->image->row0, interp); + psMetadataItem *vfItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_VARFACTOR); + if (vfItem) { + psMetadataItem *goodpixItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_GOODPIX); + psAssert(goodpixItem, "It should be where we left it!"); + psAssert(vfItem->type == PS_TYPE_F32 && goodpixItem->type == PS_TYPE_S64, + "Should be the type we said."); + + vfItem->data.F32 += varFactor * goodPixels; + goodpixItem->data.S64 += goodPixels; + } else { + psMetadataAddF32(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_VARFACTOR, 0, + "Variance factor weighted by the good pixels", varFactor * goodPixels); + psMetadataAddS64(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_GOODPIX, 0, + "Number of good pixels", goodPixels); } } - psFree(interp); - psFree(inPix); - psFree(outPix); - psFree(FP); - psFree(TP); - psFree(sky); + + if (goodPixels > 0) { + if (!pswarpTransformSources (output, input, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); + return false; + } + } + + // XXX should we not write anything out if there are no good pixels? + output->data_exists = true; + output->parent->data_exists = true; + output->parent->parent->data_exists = true; + + psLogMsg ("pswarp", 3, "warping analysis: %f sec\n", psTimerMark ("warp")); return true; Index: /trunk/pswarp/src/pswarpTransformReadout_Threaded.c =================================================================== --- /trunk/pswarp/src/pswarpTransformReadout_Threaded.c (revision 18839) +++ /trunk/pswarp/src/pswarpTransformReadout_Threaded.c (revision 18839) @@ -0,0 +1,246 @@ +# include "pswarp.h" + +# define THREADED 1 +# define SOURCE_ARRAY_BUFFER 100 // Size to grow the array of sources at a time + +// NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT +bool pswarpTransformReadout_Threaded(pmReadout *output, pmReadout *input, pmConfig *config) +{ + // XXX this implementation currently ignores the use of the region + psImage *region = NULL; + + psTimerStart ("warp"); + + // Get warp parameters + bool mdok; + int nGridX = psMetadataLookupS32(NULL, config->arguments, "GRID.NX"); + int nGridY = psMetadataLookupS32(NULL, config->arguments, "GRID.NY"); + psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.MODE"); + + // load the recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE); + psAssert (recipe, "missing recipe %s", PSWARP_RECIPE); + + // output mask bits + psMaskType maskIn = psMetadataLookupU8(&mdok, recipe, "MASK.INPUT"); + psMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config); + psMaskType maskBad = pmConfigMaskGet("BAD.WARP", config); + psAssert (mdok, "MASK.INPUT was not defined"); + + // Flux fraction for "poor" + float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC"); + + // find the output pixel range + int minX, minY, maxX, maxY; + pswarpMatchRange (&minX, &minY, &maxX, &maxY, input, output); + + // pswarpMapGridFromImage builds a set of locally-linear maps which convert the + // output coordinates to input coordinates + pswarpMapGrid *grid = pswarpMapGridFromImage (input, output, nGridX, nGridY); + + // XXX need to modify the grid based on this result and force the maxError < XXX + double maxError = pswarpMapGridMaxError (grid); + psLogMsg ("pswarp", 3, "maximum error using this grid sampling: %f\n", maxError); + + int gridX, gridY; + pswarpMapGridSetGrid (grid, minX, minY, &gridX, &gridY); +# if (0) + // XXX these asserts probably belong in the function (don't have outCol0,Row0 anyway) + assert ((int)(minX - outCol0) >= 0); + assert ((int)(maxX - outCol0) <= output->image->numCols); + assert ((int)(minY - outRow0) >= 0); + assert ((int)(maxY - outRow0) <= output->image->numRows); +# endif + + psImage *inImage = input->image; // Input image + psImage *inVar = input->weight; // Input weight map + psImage *inMask = input->mask; // Input mask + + // Interpolation options + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(interpolationMode, inImage, + inVar, inMask, maskIn, NAN, NAN, + maskBad, maskPoor, poorFrac); + + if (inVar && !output->weight) { + output->weight = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_F32); + psImageInit(output->weight, NAN); + } + if ((inMask || maskPoor || maskBad) && !output->mask) { + output->mask = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_MASK); + psImageInit(output->mask, maskBad); + } + + // total number of good pixels across all tiles (summed below) + int goodPixels = 0; + + // create jobs and supply them to the threads + for (gridY = 0; gridY < nGridX; gridY++) { + for (gridX = 0; gridX < nGridY; gridX++) { + + pswarpTransformTileArgs *args = pswarpTransformTileArgsAlloc(); + + // these items are just views to the data; they are not freed with args + args->input = input; + args->output = output; + args->grid = grid; + args->interp = interp; + args->region = region; + + args->gridX = gridX; + args->gridY = gridY; + args->goodPixels = 0; + +# if (THREADED) + // allocate a job + psThreadJob *job = psThreadJobAlloc ("PSWARP_TRANSFORM_TILE", 0); + + // construct the arguments for this job + // job is pswarpTransformTile (gridX, gridY); + psArrayAdd (job->args, 1, args); + psThreadJobAddPending (job); +# else + pswarpTransformTile (args); + goodPixels += args->goodPixels; +# endif + psFree (args); + } + } + +# if (THREADED) + // wait here for the threaded jobs to finish + if (!psThreadPoolWait ()) { + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); + return false; + } + fprintf (stderr, "success for threaded jobs\n"); + + // each job records its own goodPixel values; sum them here + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + pswarpTransformTileArgs *args = job->args->data[0]; + goodPixels += args->goodPixels; + psFree (job); + } +# endif + psFree(interp); + psFree(grid); + + + // Store the variance factor and number of good pixels + if (goodPixels > 0) { + float varFactor = psImageInterpolateVarianceFactor(inImage->numCols / 2.0 + inImage->col0, + inImage->numRows / 2.0 + inImage->row0, + interp); // Variance factor: large --> small scale + psMetadataItem *vfItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_VARFACTOR); + if (vfItem) { + psMetadataItem *goodpixItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_GOODPIX); + psAssert(goodpixItem, "It should be where we left it!"); + psAssert(vfItem->type == PS_TYPE_F32 && goodpixItem->type == PS_TYPE_S64, + "Should be the type we said."); + + vfItem->data.F32 += varFactor * goodPixels; + goodpixItem->data.S64 += goodPixels; + } else { + psMetadataAddF32(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_VARFACTOR, 0, + "Variance factor weighted by the good pixels", varFactor * goodPixels); + psMetadataAddS64(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_GOODPIX, 0, + "Number of good pixels", goodPixels); + } + } + + // Transform sources + psArray *inSources = psMetadataLookupPtr(&mdok, input->analysis, "PSPHOT.SOURCES"); // Sources in source + if (goodPixels > 0 && mdok && inSources) { + pswarpMapGrid *sourceGrid = pswarpMapGridFromImage(output, input, nGridX, nGridY); // Grid for sources + + psArray *outSources = psMemIncrRefCounter(psMetadataLookupPtr(&mdok, output->analysis, + "PSPHOT.SOURCES")); // Target sources + if (!outSources) { + outSources = psArrayAllocEmpty(SOURCE_ARRAY_BUFFER); + psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, + "Warped sources", outSources); + } + + for (int i = 0; i < inSources->n; i++) { + pmSource *source = inSources->data[i]; // Source of interest + pmModel *model = source->modelPSF; // Model for this source + float xIn, yIn; // Coordinates of source + xIn = model->params->data.F32[PM_PAR_XPOS] - inImage->col0; + yIn = model->params->data.F32[PM_PAR_YPOS] - inImage->row0; + int xGrid, yGrid; // Grid coordinates for local map + if (!pswarpMapGridSetGrid(sourceGrid, xIn + 0.5, yIn + 0.5, &xGrid, &yGrid)) { + psError(PS_ERR_UNKNOWN, false, "Unable to get grid coordinates for source at %f,%f\n", + xIn, yIn); + psFree(outSources); + psFree(sourceGrid); + return false; + } + if (xGrid < 0 || xGrid >= sourceGrid->nXpts || yGrid < 0 || yGrid >= sourceGrid->nYpts) { + // It's not even on the grid + continue; + } + + pswarpMap *map = sourceGrid->maps[xGrid][yGrid]; // Locally linear transformation + double xOut, yOut; // Output coordinates + if (!pswarpMapApply(&xOut, &yOut, map, xIn + 0.5, yIn + 0.5)) { + psError(PS_ERR_UNKNOWN, false, "Unable to transform coordinates for source at %f,%f\n", + xIn, yIn); + psFree(outSources); + psFree(sourceGrid); + return false; + } + xOut += output->image->col0 - 0.5; + yOut += output->image->row0 - 0.5; + if (xOut < minX || xOut > maxX || yOut < minY || yOut > maxY) { + // It's not in the output image + continue; + } + + // Generate the new source in the output frame + pmSource *new = pmSourceAlloc(); // New source + new->peak = pmPeakAlloc(xOut, yOut, source->peak->flux, PM_PEAK_LONE); + new->peak->flux = source->peak->flux; + new->type = PM_SOURCE_TYPE_STAR; + new->psfMag = source->psfMag; + new->errMag = source->errMag; + new->sky = source->sky; + new->skyErr = source->skyErr; + new->pixWeight = source->pixWeight; + new->modelPSF = pmModelAlloc(source->modelPSF->type); + +#if 0 + // XXX Note that this will not set the correct axes + pmPSF_AxesToModel(new->modelPSF->params->data.F32, + pmPSF_ModelToAxes(source->modelPSF->params->data.F32, 20.0)); +#endif + + // Propagate the position erorrs + float dxIn, dyIn; // Errors in input coordinates + dxIn = model->dparams->data.F32[PM_PAR_XPOS]; + dyIn = model->dparams->data.F32[PM_PAR_YPOS]; + + float dxOut, dyOut; // Errors in output coordinates + dxOut = sqrt(PS_SQR(map->Xx * dxIn) + PS_SQR(map->Xy * dyIn)); + dyOut = sqrt(PS_SQR(map->Yx * dxIn) + PS_SQR(map->Yy * dyIn)); + + new->modelPSF->params->data.F32[PM_PAR_XPOS] = xOut; + new->modelPSF->params->data.F32[PM_PAR_YPOS] = yOut; + new->modelPSF->dparams->data.F32[PM_PAR_XPOS] = dxOut; + new->modelPSF->dparams->data.F32[PM_PAR_YPOS] = dyOut; + + psArrayAdd(outSources, SOURCE_ARRAY_BUFFER, new); + psFree(new); // Drop reference + } + psFree(sourceGrid); + psFree(outSources); // Drop reference + } + + output->data_exists = true; + output->parent->data_exists = true; + output->parent->parent->data_exists = true; + + psLogMsg ("pswarp", 3, "warping analysis: %f sec\n", psTimerMark ("warp")); + + return true; +} Index: /trunk/pswarp/src/pswarpTransformReadout_Unthreaded.c =================================================================== --- /trunk/pswarp/src/pswarpTransformReadout_Unthreaded.c (revision 18839) +++ /trunk/pswarp/src/pswarpTransformReadout_Unthreaded.c (revision 18839) @@ -0,0 +1,216 @@ +# include "pswarp.h" + +#define SOURCE_ARRAY_BUFFER 100 // Size to grow the array of sources at a time + +// NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT +bool pswarpTransformReadout_Unthreaded(pmReadout *output, pmReadout *input, pmConfig *config) +{ + // XXX this implementation currently ignores the use of the region + psImage *region = NULL; + + psTimerStart ("warp"); + + // Get warp parameters + bool mdok; + int nGridX = psMetadataLookupS32(NULL, config->arguments, "GRID.NX"); + int nGridY = psMetadataLookupS32(NULL, config->arguments, "GRID.NY"); + psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.MODE"); + + // load the recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE); + psAssert (recipe, "missing recipe %s", PSWARP_RECIPE); + + // output mask bits + psMaskType maskIn = psMetadataLookupU8(&mdok, recipe, "MASK.INPUT"); + psMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config); + psMaskType maskBad = pmConfigMaskGet("BAD.WARP", config); + psAssert (mdok, "MASK.INPUT was not defined"); + + // Flux fraction for "poor" + float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC"); + + // find the output pixel range + int minX, minY, maxX, maxY; + pswarpMatchRange (&minX, &minY, &maxX, &maxY, input, output); + + // pswarpMapGridFromImage builds a set of locally-linear maps which convert the + // output coordinates to input coordinates + pswarpMapGrid *grid = pswarpMapGridFromImage (input, output, nGridX, nGridY); + + // XXX need to modify the grid based on this result and force the maxError < XXX + double maxError = pswarpMapGridMaxError (grid); + psLogMsg ("pswarp", 3, "maximum error using this grid sampling: %f\n", maxError); + + int gridX, gridY; + pswarpMapGridSetGrid (grid, minX, minY, &gridX, &gridY); +# if (0) + // XXX these asserts probably belong in the function (don't have outCol0,Row0 anyway) + assert ((int)(minX - outCol0) >= 0); + assert ((int)(maxX - outCol0) <= output->image->numCols); + assert ((int)(minY - outRow0) >= 0); + assert ((int)(maxY - outRow0) <= output->image->numRows); +# endif + + psImage *inImage = input->image; // Input image + psImage *inVar = input->weight; // Input weight map + psImage *inMask = input->mask; // Input mask + + // Interpolation options + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(interpolationMode, inImage, + inVar, inMask, maskIn, NAN, NAN, + maskBad, maskPoor, poorFrac); + + if (inVar && !output->weight) { + output->weight = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_F32); + psImageInit(output->weight, NAN); + } + if ((inMask || maskPoor || maskBad) && !output->mask) { + output->mask = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_MASK); + psImageInit(output->mask, maskBad); + } + + // total number of good pixels across all tiles (summed below) + int goodPixels = 0; + + // (almost) the same structure as the threaded version + for (gridY = 0; gridY < nGridX; gridY++) { + for (gridX = 0; gridX < nGridY; gridX++) { + + pswarpTransformTileArgs *args = pswarpTransformTileArgsAlloc(); + + // these items are just views to the data; they are not freed with args + args->input = input; + args->output = output; + args->grid = grid; + args->interp = interp; + args->region = region; + + args->gridX = gridX; + args->gridY = gridY; + args->goodPixels = 0; + + pswarpTransformTile (args); + goodPixels += args->goodPixels; + psFree (args); + } + } + psFree(interp); + psFree(grid); + + // Store the variance factor and number of good pixels + if (goodPixels > 0) { + float varFactor = psImageInterpolateVarianceFactor(inImage->numCols / 2.0 + inImage->col0, + inImage->numRows / 2.0 + inImage->row0, + interp); // Variance factor: large --> small scale + psMetadataItem *vfItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_VARFACTOR); + if (vfItem) { + psMetadataItem *goodpixItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_GOODPIX); + psAssert(goodpixItem, "It should be where we left it!"); + psAssert(vfItem->type == PS_TYPE_F32 && goodpixItem->type == PS_TYPE_S64, + "Should be the type we said."); + + vfItem->data.F32 += varFactor * goodPixels; + goodpixItem->data.S64 += goodPixels; + } else { + psMetadataAddF32(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_VARFACTOR, 0, + "Variance factor weighted by the good pixels", varFactor * goodPixels); + psMetadataAddS64(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_GOODPIX, 0, + "Number of good pixels", goodPixels); + } + } + + // Transform sources + psArray *inSources = psMetadataLookupPtr(&mdok, input->analysis, "PSPHOT.SOURCES"); // Sources in source + if (goodPixels > 0 && mdok && inSources) { + pswarpMapGrid *sourceGrid = pswarpMapGridFromImage(output, input, nGridX, nGridY); // Grid for sources + + psArray *outSources = psMemIncrRefCounter(psMetadataLookupPtr(&mdok, output->analysis, + "PSPHOT.SOURCES")); // Target sources + if (!outSources) { + outSources = psArrayAllocEmpty(SOURCE_ARRAY_BUFFER); + psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, + "Warped sources", outSources); + } + + for (int i = 0; i < inSources->n; i++) { + pmSource *source = inSources->data[i]; // Source of interest + pmModel *model = source->modelPSF; // Model for this source + float xIn, yIn; // Coordinates of source + xIn = model->params->data.F32[PM_PAR_XPOS] - inImage->col0; + yIn = model->params->data.F32[PM_PAR_YPOS] - inImage->row0; + int xGrid, yGrid; // Grid coordinates for local map + if (!pswarpMapGridSetGrid(sourceGrid, xIn + 0.5, yIn + 0.5, &xGrid, &yGrid)) { + psError(PS_ERR_UNKNOWN, false, "Unable to get grid coordinates for source at %f,%f\n", + xIn, yIn); + psFree(outSources); + psFree(sourceGrid); + return false; + } + if (xGrid < 0 || xGrid >= sourceGrid->nXpts || yGrid < 0 || yGrid >= sourceGrid->nYpts) { + // It's not even on the grid + continue; + } + + pswarpMap *map = sourceGrid->maps[xGrid][yGrid]; // Locally linear transformation + double xOut, yOut; // Output coordinates + if (!pswarpMapApply(&xOut, &yOut, map, xIn + 0.5, yIn + 0.5)) { + psError(PS_ERR_UNKNOWN, false, "Unable to transform coordinates for source at %f,%f\n", + xIn, yIn); + psFree(outSources); + psFree(sourceGrid); + return false; + } + xOut += output->image->col0 - 0.5; + yOut += output->image->row0 - 0.5; + if (xOut < minX || xOut > maxX || yOut < minY || yOut > maxY) { + // It's not in the output image + continue; + } + + // Generate the new source in the output frame + pmSource *new = pmSourceAlloc(); // New source + new->peak = pmPeakAlloc(xOut, yOut, source->peak->flux, PM_PEAK_LONE); + new->peak->flux = source->peak->flux; + new->type = PM_SOURCE_TYPE_STAR; + new->psfMag = source->psfMag; + new->errMag = source->errMag; + new->sky = source->sky; + new->skyErr = source->skyErr; + new->pixWeight = source->pixWeight; + new->modelPSF = pmModelAlloc(source->modelPSF->type); + +#if 0 + // XXX Note that this will not set the correct axes + pmPSF_AxesToModel(new->modelPSF->params->data.F32, + pmPSF_ModelToAxes(source->modelPSF->params->data.F32, 20.0)); +#endif + + // Propagate the position erorrs + float dxIn, dyIn; // Errors in input coordinates + dxIn = model->dparams->data.F32[PM_PAR_XPOS]; + dyIn = model->dparams->data.F32[PM_PAR_YPOS]; + + float dxOut, dyOut; // Errors in output coordinates + dxOut = sqrt(PS_SQR(map->Xx * dxIn) + PS_SQR(map->Xy * dyIn)); + dyOut = sqrt(PS_SQR(map->Yx * dxIn) + PS_SQR(map->Yy * dyIn)); + + new->modelPSF->params->data.F32[PM_PAR_XPOS] = xOut; + new->modelPSF->params->data.F32[PM_PAR_YPOS] = yOut; + new->modelPSF->dparams->data.F32[PM_PAR_XPOS] = dxOut; + new->modelPSF->dparams->data.F32[PM_PAR_YPOS] = dyOut; + + psArrayAdd(outSources, SOURCE_ARRAY_BUFFER, new); + psFree(new); // Drop reference + } + psFree(sourceGrid); + psFree(outSources); // Drop reference + } + + output->data_exists = true; + output->parent->data_exists = true; + output->parent->parent->data_exists = true; + + psLogMsg ("pswarp", 3, "warping analysis: %f sec\n", psTimerMark ("warp")); + + return true; +} Index: /trunk/pswarp/src/pswarpTransformSources.c =================================================================== --- /trunk/pswarp/src/pswarpTransformSources.c (revision 18839) +++ /trunk/pswarp/src/pswarpTransformSources.c (revision 18839) @@ -0,0 +1,110 @@ +# include "pswarp.h" + +# define SOURCE_ARRAY_BUFFER 100 // Size to grow the array of sources at a time + +// NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT +bool pswarpTransformSources(pmReadout *output, pmReadout *input, pmConfig *config) +{ + + // find the output pixel range + int minX, minY, maxX, maxY; + pswarpMatchRange (&minX, &minY, &maxX, &maxY, input, output); + + // Get warp parameters + bool mdok; + int nGridX = psMetadataLookupS32(NULL, config->arguments, "GRID.NX"); + int nGridY = psMetadataLookupS32(NULL, config->arguments, "GRID.NY"); + + // Transform sources + psArray *inSources = psMetadataLookupPtr(&mdok, input->analysis, "PSPHOT.SOURCES"); // Sources in source + if (!inSources) return true; + + pswarpMapGrid *sourceGrid = pswarpMapGridFromImage(output, input, nGridX, nGridY); // Grid for sources + + psArray *outSources = psMemIncrRefCounter(psMetadataLookupPtr(&mdok, output->analysis, "PSPHOT.SOURCES")); // Target sources + if (!outSources) { + outSources = psArrayAllocEmpty(SOURCE_ARRAY_BUFFER); + psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, + "Warped sources", outSources); + } + + // XXX it is probably not necessary to use the locally linear transformations we should + // be using the full astrometry: since there are 100 - 1000 fewer sources than pixels, + // this does not cost us so much time. + for (int i = 0; i < inSources->n; i++) { + pmSource *source = inSources->data[i]; // Source of interest + pmModel *model = source->modelPSF; // Model for this source + float xIn, yIn; // Coordinates of source + xIn = model->params->data.F32[PM_PAR_XPOS] - input->image->col0; + yIn = model->params->data.F32[PM_PAR_YPOS] - input->image->row0; + + int xGrid, yGrid; // Grid coordinates for local map + if (!pswarpMapGridSetGrid(sourceGrid, xIn + 0.5, yIn + 0.5, &xGrid, &yGrid)) { + psError(PS_ERR_UNKNOWN, false, "Unable to get grid coordinates for source at %f,%f\n", + xIn, yIn); + psFree(outSources); + psFree(sourceGrid); + return false; + } + if (xGrid < 0 || xGrid >= sourceGrid->nXpts || yGrid < 0 || yGrid >= sourceGrid->nYpts) { + // It's not even on the grid + // XXX how can this happen? + continue; + } + + pswarpMap *map = sourceGrid->maps[xGrid][yGrid]; // Locally linear transformation + double xOut, yOut; // Output coordinates + if (!pswarpMapApply(&xOut, &yOut, map, xIn + 0.5, yIn + 0.5)) { + psError(PS_ERR_UNKNOWN, false, "Unable to transform coordinates for source at %f,%f\n", + xIn, yIn); + psFree(outSources); + psFree(sourceGrid); + return false; + } + xOut += output->image->col0 - 0.5; + yOut += output->image->row0 - 0.5; + if (xOut < minX || xOut > maxX || yOut < minY || yOut > maxY) { + // It's not in the output image + continue; + } + + // Generate the new source in the output frame + pmSource *new = pmSourceAlloc(); // New source + new->peak = pmPeakAlloc(xOut, yOut, source->peak->flux, PM_PEAK_LONE); + new->peak->flux = source->peak->flux; + new->type = PM_SOURCE_TYPE_STAR; + new->psfMag = source->psfMag; + new->errMag = source->errMag; + new->sky = source->sky; + new->skyErr = source->skyErr; + new->pixWeight = source->pixWeight; + new->modelPSF = pmModelAlloc(source->modelPSF->type); + +#if 0 + // XXX Note that this will not set the correct axes + pmPSF_AxesToModel(new->modelPSF->params->data.F32, + pmPSF_ModelToAxes(source->modelPSF->params->data.F32, 20.0)); +#endif + + // Propagate the position erorrs + float dxIn, dyIn; // Errors in input coordinates + dxIn = model->dparams->data.F32[PM_PAR_XPOS]; + dyIn = model->dparams->data.F32[PM_PAR_YPOS]; + + float dxOut, dyOut; // Errors in output coordinates + dxOut = sqrt(PS_SQR(map->Xx * dxIn) + PS_SQR(map->Xy * dyIn)); + dyOut = sqrt(PS_SQR(map->Yx * dxIn) + PS_SQR(map->Yy * dyIn)); + + new->modelPSF->params->data.F32[PM_PAR_XPOS] = xOut; + new->modelPSF->params->data.F32[PM_PAR_YPOS] = yOut; + new->modelPSF->dparams->data.F32[PM_PAR_XPOS] = dxOut; + new->modelPSF->dparams->data.F32[PM_PAR_YPOS] = dyOut; + + psArrayAdd(outSources, SOURCE_ARRAY_BUFFER, new); + psFree(new); // Drop reference + } + psFree(sourceGrid); + psFree(outSources); // Drop reference + return true; +} + Index: /trunk/pswarp/src/pswarpTransformTile.c =================================================================== --- /trunk/pswarp/src/pswarpTransformTile.c (revision 18839) +++ /trunk/pswarp/src/pswarpTransformTile.c (revision 18839) @@ -0,0 +1,103 @@ +# include "pswarp.h" + +void pswarpTransformTileArgsFree (pswarpTransformTileArgs *args) { + if (!args) return; + return; +} + +pswarpTransformTileArgs *pswarpTransformTileArgsAlloc() { + + pswarpTransformTileArgs *args = (pswarpTransformTileArgs *)psAlloc(sizeof(pswarpTransformTileArgs)); + psMemSetDeallocator(args, (psFreeFunc)pswarpTransformTileArgsFree); + + args->input = NULL; + args->output = NULL; + args->grid = NULL; + args->interp = NULL; + args->region = NULL; + + args->gridX = 0; + args->gridY = 0; + + args->goodPixels = 0; + + return args; +} + +bool pswarpTransformTile (pswarpTransformTileArgs *args) { + + // int inCol0 = args->input->image->col0; + // int inRow0 = args->input->image->row0; + int inNcol = args->input->image->numCols; + int inNrow = args->input->image->numRows; + + int outCol0 = args->output->image->col0; + int outRow0 = args->output->image->row0; + // int outNcol = args->output->image->numCols; + // int outNrow = args->output->image->numRows; + + // get the coordinate range for this grid tile + psPlane minPt, maxPt; + pswarpMapGridCoordRange (args->grid, args->gridX, args->gridY, &minPt, &maxPt); + + psF32 **outImageData = (args->output->image) ? args->output->image->data.F32 : NULL; + psF32 **outVarData = (args->output->weight) ? args->output->weight->data.F32 : NULL; + psMaskType **outMaskData = (args->output->mask) ? args->output->mask->data.U8 : NULL; + + pswarpMap *map = args->grid->maps[args->gridX][args->gridY]; + + psImage *region = args->region; + + double xInRaw, yInRaw; + + // output values for this pixel + double imageValue; + double varValue; + psMaskType maskValue; + + // Iterate over the output image pixels (parent frame) + long goodPixels = 0; // Number of input pixels landing on the output image + for (int y = minPt.y; y < maxPt.y; y++) { + for (int x = minPt.x; x < maxPt.x; x++) { + + // Only transform those pixels requested + if (region && region->data.U8[y][x]) continue; + + // pswarpMapApply converts the output coordinate (x,y) to the input coordinate. + // both are in the parent frames of the input and output images. + pswarpMapApply (&xInRaw, &yInRaw, map, x + 0.5, y + 0.5); + + double xIn = xInRaw - outCol0; // Position on input image + double yIn = yInRaw - outRow0; // Position on input image + + if (xIn < 0) continue; + if (yIn >= inNcol) continue; + if (yIn < 0) continue; + if (yIn >= inNrow) continue; + + goodPixels++; + + // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates + if (!psImageInterpolate(&imageValue, &varValue, &maskValue, xIn, yIn, args->interp)) { + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); + return false; + } + + int xOut = x - outCol0; // Position on output image + int yOut = y - outRow0; // Position on output image + + // not all images need be transformed + if (outImageData) { + outImageData[yOut][xOut] = imageValue; + } + if (outVarData) { + outVarData[yOut][xOut] = varValue; + } + if (outMaskData) { + outMaskData[yOut][xOut] = maskValue; + } + } + } + args->goodPixels = goodPixels; + return true; +}