Index: /branches/eam_branch_20080719/pswarp/src/Makefile.am =================================================================== --- /branches/eam_branch_20080719/pswarp/src/Makefile.am (revision 18751) +++ /branches/eam_branch_20080719/pswarp/src/Makefile.am (revision 18752) @@ -16,6 +16,9 @@ pswarpPixelFraction.c \ pswarpSetMaskBits.c \ + pswarpThreadLauncher.c \ pswarpTransformReadout_Opt.c \ - pswarpThreadLauncher.c \ + pswarpTransformReadout_Threaded.c \ + pswarpTransformReadout_Unthreaded.c \ + pswarpTransformTile.c \ pswarpVersion.c Index: /branches/eam_branch_20080719/pswarp/src/pswarp.c =================================================================== --- /branches/eam_branch_20080719/pswarp/src/pswarp.c (revision 18751) +++ /branches/eam_branch_20080719/pswarp/src/pswarp.c (revision 18752) @@ -18,30 +18,7 @@ // create the thread pool with 4 threads, supplying our thread launcher function + // XXX need to determine the number of threads from the config data psThreadPoolInit (4, &pswarpThreadLauncher); - srand48(0); - - for (int i = 0; i < 10; i++) { - // allocate a job which takes two arguments - // XXX probably don't need to pre-specify the number of args - psThreadJob *job = psThreadJobAlloc ("THREAD_TEST", 2); - - // construct the arguments for this job - char *arg1 = NULL; - char *arg2 = NULL; - psStringAppend (&arg1, "arg1.%d", i); - psStringAppend (&arg2, "arg2.%d", i); - psArrayAdd (job->args, 1, arg1); - psArrayAdd (job->args, 1, arg2); - - psThreadJobAddToQueue (job); - } - if (!psThreadPoolWait ()) { - fprintf (stderr, "error in a threaded job\n"); - exit (1); - } - fprintf (stderr, "success for threaded jobs\n"); - exit (0); - // model inits are needed in pmSourceIO // models defined in psphot/src/models are not available in psastro Index: /branches/eam_branch_20080719/pswarp/src/pswarp.h =================================================================== --- /branches/eam_branch_20080719/pswarp/src/pswarp.h (revision 18751) +++ /branches/eam_branch_20080719/pswarp/src/pswarp.h (revision 18752) @@ -36,4 +36,26 @@ } 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; + +bool pswarpTransformReadout_Unthreaded(pmReadout *output, pmReadout *input, pmConfig *config); +bool pswarpTransformReadout_Threaded(pmReadout *output, pmReadout *input, pmConfig *config); + +pswarpTransformTileArgs *pswarpTransformTileArgsAlloc(); +bool pswarpTransformTile (pswarpTransformTileArgs *args); + pmConfig *pswarpArguments (int argc, char **argv); bool pswarpOptions(pmConfig *config); Index: /branches/eam_branch_20080719/pswarp/src/pswarpThreadLauncher.c =================================================================== --- /branches/eam_branch_20080719/pswarp/src/pswarpThreadLauncher.c (revision 18751) +++ /branches/eam_branch_20080719/pswarp/src/pswarpThreadLauncher.c (revision 18752) @@ -33,5 +33,5 @@ // we have to lock here so the job queue cannot be empty yet no threads busy psThreadLock(); - while ((job = psThreadJobGet ()) == NULL) { + while ((job = psThreadJobGetPending ()) == NULL) { psThreadUnlock(); usleep (10000); @@ -50,4 +50,16 @@ continue; } + + // 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: /branches/eam_branch_20080719/pswarp/src/pswarpTransformReadout_Threaded.c =================================================================== --- /branches/eam_branch_20080719/pswarp/src/pswarpTransformReadout_Threaded.c (revision 18751) +++ /branches/eam_branch_20080719/pswarp/src/pswarpTransformReadout_Threaded.c (revision 18752) @@ -15,12 +15,9 @@ int nGridX = psMetadataLookupS32(NULL, config->arguments, "GRID.NX"); int nGridY = psMetadataLookupS32(NULL, config->arguments, "GRID.NY"); - psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments, - "INTERPOLATION.MODE"); + psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.MODE"); + // load the recipe psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE); - if (!recipe) { - psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE); - return false; - } + psAssert (recipe, "missing recipe %s", PSWARP_RECIPE); // output mask bits @@ -30,15 +27,6 @@ psAssert (mdok, "MASK.INPUT was not defined"); - float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC"); // Flux fraction for "poor" - - // we need to apply the offset to convert parent coordinates to child coordinates for - // psImagePixelInterpolate below - int inCol0 = input->image->col0; - int inRow0 = input->image->row0; - int outCol0 = output->image->col0; - int outRow0 = output->image->row0; - - // we might want to do the rectangular regions outside of the selection independently - // psImageInit (output->image, NAN); + // Flux fraction for "poor" + float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC"); // find the output pixel range @@ -54,16 +42,13 @@ psLogMsg ("pswarp", 3, "maximum error using this grid sampling: %f\n", maxError); - int gridX, gridY, nextGridX, nextGridY; + int gridX, gridY; pswarpMapGridSetGrid (grid, minX, minY, &gridX, &gridY); - nextGridY = pswarpMapGridNextGrid_Y (grid, gridY); - pswarpMap *map = grid->maps[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); - - int gridXo = gridX; - int nextGridXo = pswarpMapGridNextGrid_X (grid, gridX); +# endif psImage *inImage = input->image; // Input image @@ -76,101 +61,65 @@ maskBad, maskPoor, poorFrac); - psPlane *inPix = psPlaneAlloc(); // Coordinates on the input detector - psF32 **outImageData = output->image->data.F32; // Output image pixels - psF32 **outVarData = NULL; // Output variance pixels - psMaskType **outMaskData = NULL; // Output mask pixels - if (inVar) { - if (!output->weight) { - output->weight = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_F32); - psImageInit(output->weight, NAN); - } - outVarData = output->weight->data.F32; - } - if (inMask || maskPoor || maskBad) { - if (!output->mask) { - output->mask = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_MASK); - psImageInit(output->mask, maskBad); - } - outMaskData = output->mask->data.PS_TYPE_MASK_DATA; - } - - // create N threads - for (gridY = gridYmin; gridY < gridYmax; gridY++) { - for (gridX = gridXmin; gridX < gridXmax; gridX++) { - pswarpTransformTile (gridX, gridY); - // send function and arguments to thread manager - // threadFunction (pswarpTransformTile) - // while (threadFull()) { - // usleep (50000); - // } + 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; + + // 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); + psFree (args); } } - - - // Iterate over the output image pixels (parent frame) - long goodPixels = 0; // Number of input pixels landing on the output image - for (int y = minY; y < maxY; y++) { - if (y >= nextGridY) { - gridY ++; - nextGridY = pswarpMapGridNextGrid_Y (grid, gridY); - map = grid->maps[gridX][gridY]; - } - - gridX = gridXo; - nextGridX = nextGridXo; - map = grid->maps[gridX][gridY]; - int yOut = y - outRow0; // Position on image - for (int x = minX; x < maxX; x++) { - if (x >= nextGridX) { - gridX ++; - nextGridX = pswarpMapGridNextGrid_X (grid, gridX); - map = grid->maps[gridX][gridY]; - } - - // 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 (&inPix->x, &inPix->y, map, x + 0.5, y + 0.5); - - if (inPix->x - inCol0 < 0) continue; - if (inPix->x - inCol0 >= inImage->numCols) continue; - if (inPix->y - inRow0 < 0) continue; - if (inPix->y - inRow0 >= inImage->numRows) continue; - - goodPixels++; - - // XXX include mask - // XXX apply scale and offset? - // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates - double imageValue, varValue; - psMaskType maskValue; - if (!psImageInterpolate(&imageValue, &varValue, &maskValue, - inPix->x - inCol0, inPix->y - inRow0, interp)) { - psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); - psFree(interp); - psFree(inPix); - psFree(grid); - return false; - } - int xOut = x - outCol0; // Position on image - outImageData[yOut][xOut] = imageValue; - if (inVar) { - outVarData[yOut][xOut] = varValue; - } - if (outMaskData) { - outMaskData[yOut][xOut] = maskValue; - } - } + // 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); } psFree(interp); - psFree(inPix); psFree(grid); // Store the variance factor and number of good pixels if (goodPixels > 0) { - float varFactor = psImageInterpolateVarianceFactor(inImage->numCols / 2.0 + inCol0, - inImage->numRows / 2.0 + inRow0, + 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); @@ -208,6 +157,6 @@ pmModel *model = source->modelPSF; // Model for this source float xIn, yIn; // Coordinates of source - xIn = model->params->data.F32[PM_PAR_XPOS] - inCol0; - yIn = model->params->data.F32[PM_PAR_YPOS] - inRow0; + 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)) { @@ -232,6 +181,6 @@ return false; } - xOut += outCol0 - 0.5; - yOut += outRow0 - 0.5; + 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 Index: /branches/eam_branch_20080719/pswarp/src/pswarpTransformReadout_Unthreaded.c =================================================================== --- /branches/eam_branch_20080719/pswarp/src/pswarpTransformReadout_Unthreaded.c (revision 18752) +++ /branches/eam_branch_20080719/pswarp/src/pswarpTransformReadout_Unthreaded.c (revision 18752) @@ -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: /branches/eam_branch_20080719/pswarp/src/pswarpTransformTile.c =================================================================== --- /branches/eam_branch_20080719/pswarp/src/pswarpTransformTile.c (revision 18751) +++ /branches/eam_branch_20080719/pswarp/src/pswarpTransformTile.c (revision 18752) @@ -1,26 +1,64 @@ +# include "pswarp.h" -// XXX need to determine all of the needed inputs and put them in a structure +void pswarpTransformTileArgsFree (pswarpTransformTileArgs *args) { + if (!args) return; + return; +} -typedef struct { - int gridX; - int gridY; -} pswarpTransformTileOpts; +pswarpTransformTileArgs *pswarpTransformTileArgsAlloc() { -pswarpTransformTile (int gridX, int gridY) { + pswarpTransformTileArgs *args = (pswarpTransformTileArgs *)psAlloc(sizeof(pswarpTransformTileArgs)); + psMemSetDeallocator(args, (psFreeFunc)pswarpTransformTileArgsFree); - minX = pswarpMapGridNextGrid_Y (grid, gridY); - minY = pswarpMapGridNextGrid_X (grid, gridX); + args->input = NULL; + args->output = NULL; + args->grid = NULL; + args->interp = NULL; + args->region = NULL; - maxX = pswarpMapGridNextGrid_Y (grid, gridY); - maxY = pswarpMapGridNextGrid_X (grid, gridX); + args->gridX = 0; + args->gridY = 0; - map = grid->maps[gridX][gridY]; + args->goodPixels = 0; - psPlane *inPix = psPlaneAlloc(); // Coordinates on the input detector + 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; + + // XXX these are wrong (min or max is wrong) + int minX = pswarpMapGridNextGrid_Y (args->grid, args->gridY); + int minY = pswarpMapGridNextGrid_X (args->grid, args->gridX); + int maxX = pswarpMapGridNextGrid_Y (args->grid, args->gridY); + int maxY = pswarpMapGridNextGrid_X (args->grid, args->gridX); + + psF32 **outImageData = args->output->image->data.F32; + psF32 **outVarData = args->output->weight->data.F32; + psMaskType **outMaskData = args->output->mask->data.U8; + + 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 = minY; y < maxY; y++) { - int yOut = y - outRow0; // Position on image for (int x = minX; x < maxX; x++) { @@ -30,29 +68,30 @@ // pswarpMapApply converts the output coordinate (x,y) to the input coordinate. // both are in the parent frames of the input and output images. - pswarpMapApply (&inPix->x, &inPix->y, map, x + 0.5, y + 0.5); + pswarpMapApply (&xInRaw, &yInRaw, map, x + 0.5, y + 0.5); - if (inPix->x - inCol0 < 0) continue; - if (inPix->x - inCol0 >= inImage->numCols) continue; - if (inPix->y - inRow0 < 0) continue; - if (inPix->y - inRow0 >= inImage->numRows) continue; + 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++; - // XXX include mask - // XXX apply scale and offset? // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates - double imageValue, varValue; - psMaskType maskValue; - if (!psImageInterpolate(&imageValue, &varValue, &maskValue, - inPix->x - inCol0, inPix->y - inRow0, interp)) { + if (!psImageInterpolate(&imageValue, &varValue, &maskValue, xIn, yIn, args->interp)) { psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); - psFree(interp); - psFree(inPix); - psFree(grid); return false; } - int xOut = x - outCol0; // Position on image - outImageData[yOut][xOut] = imageValue; - if (inVar) { + + 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; } @@ -62,5 +101,5 @@ } } - - psFree(inPix); + args->goodPixels = goodPixels; + return true; }