Index: trunk/pswarp/src/pswarpTransformReadout.c =================================================================== --- trunk/pswarp/src/pswarpTransformReadout.c (revision 12744) +++ 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;