Index: /trunk/pswarp/src/pswarpTransformReadout.c =================================================================== --- /trunk/pswarp/src/pswarpTransformReadout.c (revision 12743) +++ /trunk/pswarp/src/pswarpTransformReadout.c (revision 12744) @@ -3,5 +3,5 @@ bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config) { - // XXX this implementation currently ignores the use of the region + // XXX this implementation currently ignores the use of the region psImage *region = NULL; pmCell *cell = NULL; @@ -13,10 +13,10 @@ int outNy = output->image->numRows; - psPlane *inPix = psPlaneAlloc(); // Coordinates on the input detector - psPlane *outPix = psPlaneAlloc(); // Coordinates on the output detector + psPlane *inPix = psPlaneAlloc(); // Coordinates on the input detector + psPlane *outPix = psPlaneAlloc(); // Coordinates on the output detector - psPlane *FP = psPlaneAlloc(); // Coordinates on the focal plane - psPlane *TP = psPlaneAlloc(); // Coordinates on the tangent plane - psSphere *sky = psSphereAlloc(); // Coordinates on the sky + psPlane *FP = psPlaneAlloc(); // Coordinates on the focal plane + psPlane *TP = psPlaneAlloc(); // Coordinates on the tangent plane + psSphere *sky = psSphereAlloc(); // Coordinates on the sky cell = input->parent; @@ -31,48 +31,49 @@ psImage *inImage = input->image; + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR, inImage, + NULL, NULL, 0, NAN, NAN, 0, 0, 0.0); + // 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; + for (int x = 0; x < outNx; x++) { + // Only transform those pixels requested + if (region && region->data.U8[y][x]) continue; - // XXX double check this 1/2 pixel offset - outPix->x = (double)x + 0.5; - outPix->y = (double)y + 0.5; + // XXX double check this 1/2 pixel offset + outPix->x = (double)x + 0.5; + outPix->y = (double)y + 0.5; - psPlaneTransformApply(FP, chipOutput->toFPA, outPix); - psPlaneTransformApply (TP, fpaOutput->toTPA, FP); - psDeproject (sky, TP, fpaOutput->toSky); - - psProject (TP, sky, fpaInput->toSky); - psPlaneTransformApply (FP, fpaInput->fromTPA, TP); - psPlaneTransformApply (inPix, chipInput->fromFPA, FP); + psPlaneTransformApply(FP, chipOutput->toFPA, outPix); + psPlaneTransformApply (TP, fpaOutput->toTPA, FP); + psDeproject (sky, TP, fpaOutput->toSky); - // XXX get interpolation method from the recipe - outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, NULL, 1, NAN, PS_INTERPOLATE_BILINEAR); - // outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, mask, 1, NAN, PS_INTERPOLATE_BILINEAR); - // modify zero and scale? - } + psProject (TP, sky, fpaInput->toSky); + psPlaneTransformApply (FP, fpaInput->fromTPA, TP); + psPlaneTransformApply (inPix, chipInput->fromFPA, FP); + + // 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; + } + + outData[y][x] = value; + // modify zero and scale? + } } - - psFree (inPix); - psFree (outPix); - psFree (FP); - psFree (TP); - psFree (sky); + psFree(interp); + psFree(inPix); + psFree(outPix); + psFree(FP); + psFree(TP); + psFree(sky); return true; } - -# if (0) - if (error) { - // Error is actually the variance - outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error, - detector->x, - detector->y, - mask, 1, NAN); - } - if (error) { - outError->data.F32[y][x] = outError->data.F32[y][x] / SQUARE(scale); - } -# endif Index: /trunk/pswarp/src/pswarpTransformReadout_Opt.c =================================================================== --- /trunk/pswarp/src/pswarpTransformReadout_Opt.c (revision 12743) +++ /trunk/pswarp/src/pswarpTransformReadout_Opt.c (revision 12744) @@ -1,4 +1,3 @@ # include "pswarp.h" -psImageInterpolateMode psImageInterpolateModeFromString (char *name); // NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT @@ -11,5 +10,5 @@ psImageInterpolateMode interpolationMode = PS_INTERPOLATE_NONE; - // XXX this implementation currently ignores the use of the region + // XXX this implementation currently ignores the use of the region psImage *region = NULL; @@ -19,8 +18,8 @@ psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE); if (!recipe) { - psError(PSWARP_ERR_CONFIG, true, "Can't find PSWARP recipe!\n"); - return false; + psError(PSWARP_ERR_CONFIG, true, "Can't find PSWARP recipe!\n"); + return false; } - + int nGridX = psMetadataLookupS32 (&status, recipe, "GRID.NX"); if (!status) nGridX = 128; @@ -30,15 +29,15 @@ char *name = psMetadataLookupStr (&status, recipe, "INTERPOLATION.MODE"); if (!name) { - interpolationMode = PS_INTERPOLATE_BILINEAR; - psLogMsg ("pswarp", 3, "defaulting to bilinear interpolation\n"); + interpolationMode = PS_INTERPOLATE_BILINEAR; + psLogMsg ("pswarp", 3, "defaulting to bilinear interpolation\n"); } else { - interpolationMode = psImageInterpolateModeFromString (name); - if (interpolationMode == PS_INTERPOLATE_NONE) { - interpolationMode = PS_INTERPOLATE_BILINEAR; - psLogMsg ("pswarp", 3, "unknown interpolation mode %s, defaulting to bilinear interpolation\n", name); - } + interpolationMode = psImageInterpolateModeFromString (name); + if (interpolationMode == PS_INTERPOLATE_NONE) { + interpolationMode = PS_INTERPOLATE_BILINEAR; + psLogMsg ("pswarp", 3, "unknown interpolation mode %s, defaulting to bilinear interpolation\n", name); + } } - psPlane *inPix = psPlaneAlloc(); // Coordinates on the input detector + psPlane *inPix = psPlaneAlloc(); // Coordinates on the input detector psImage *inImage = input->image; psF32 **outData = output->image->data.F32; @@ -65,5 +64,5 @@ psLogMsg ("pswarp", 3, "maximum error using this grid sampling: %f\n", maxError); - pswarpMapGridSetGrid (grid, minX, minY, &gridX, &gridY); + pswarpMapGridSetGrid (grid, minX, minY, &gridX, &gridY); nextGridY = pswarpMapGridNextGrid_Y (grid, gridY); map = grid->maps[gridX][gridY]; @@ -77,41 +76,53 @@ nextGridXo = pswarpMapGridNextGrid_X (grid, gridX); + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR, inImage, + NULL, NULL, 0, NAN, NAN, 0, 0, 0.0); + // Iterate over the output image pixels (parent frame) for (int y = minY; y < maxY; y++) { - if (y >= nextGridY) { - gridY ++; - nextGridY = pswarpMapGridNextGrid_Y (grid, gridY); - map = grid->maps[gridX][gridY]; - } + if (y >= nextGridY) { + gridY ++; + nextGridY = pswarpMapGridNextGrid_Y (grid, gridY); + map = grid->maps[gridX][gridY]; + } - gridX = gridXo; - nextGridX = nextGridXo; - map = grid->maps[gridX][gridY]; - for (int x = minX; x < maxX; x++) { - if (x >= nextGridX) { - gridX ++; - nextGridX = pswarpMapGridNextGrid_X (grid, gridX); - map = grid->maps[gridX][gridY]; - } + gridX = gridXo; + nextGridX = nextGridXo; + map = grid->maps[gridX][gridY]; + 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; + // 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); + // 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; + 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; - // XXX include mask - // XXX apply scale and offset? - // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates - outData[y-outRow0][x-outCol0] = (psF32)psImagePixelInterpolate(inImage, inPix->x - inCol0, inPix->y - inRow0, NULL, 1, NAN, interpolationMode); - } + // XXX include mask + // XXX apply scale and offset? + // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates + double value; + if (!psImageInterpolate(&value, NULL, NULL, inPix->x - inCol0, inPix->y - inRow0, interp)) { + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); + psFree(interp); + psFree (inPix); + psFree (grid); + return false; + } + outData[y-outRow0][x-outCol0] = value; + } } + psFree(interp); psFree (inPix); psFree (grid); Index: /trunk/stac/src/stacRejection.c =================================================================== --- /trunk/stac/src/stacRejection.c (revision 12743) +++ /trunk/stac/src/stacRejection.c (revision 12744) @@ -111,4 +111,8 @@ #endif + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR, + reject, NULL, NULL, 0, 0.0, 0.0, + 0, 0, 0.0); + // Transform the mask // Optimisation option is to only transform the pixels that have been rejected in the output, @@ -125,6 +129,14 @@ inCoords->y = (double)y + 0.5; (void)psPlaneTransformApply(outCoords, map, inCoords); - float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y, - NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR); + double maskVal; + if (!psImageInterpolate(&maskVal, NULL, NULL, outCoords->x, outCoords->y, interp)) { + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); + psFree(grads); + psFree(gradsMask); + psFree(inCoords); + psFree(outCoords); + return NULL; + } + #ifdef TESTING rejmap->data.F32[y][x] = maskVal; Index: /trunk/stac/src/stacTransform.c =================================================================== --- /trunk/stac/src/stacTransform.c (revision 12743) +++ /trunk/stac/src/stacTransform.c (revision 12744) @@ -193,4 +193,8 @@ } + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR, + image, error, mask, 1, NAN, NAN, + 1, 0, 0.0); + // Iterate over the output image pixels for (int y = 0; y < outny; y++) { @@ -204,13 +208,18 @@ // Change PS_INTERPOLATE_BILINEAR to best available technique. - outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, detector->x, - detector->y, mask, 1, NAN, - PS_INTERPOLATE_BILINEAR); + double imageValue, varianceValue; + if (!psImageInterpolate(&imageValue, &varianceValue, NULL, + detector->x, detector->y, interp)) { + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); + psFree(interp); + psFree(detector); + psFree(sky); + return false; + } + outImage->data.F32[y][x] = imageValue; + if (error) { // Error is actually the variance - outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error, - detector->x, - detector->y, - mask, 1, NAN); + outError->data.F32[y][x] = varianceValue; } @@ -224,4 +233,5 @@ } } // Iterating over output pixels + psFree(interp); } // Iterating over images