Index: /branches/czw_branch/20120906/psLib/src/imageops/psImageInterpolate.c =================================================================== --- /branches/czw_branch/20120906/psLib/src/imageops/psImageInterpolate.c (revision 34469) +++ /branches/czw_branch/20120906/psLib/src/imageops/psImageInterpolate.c (revision 34470) @@ -947,31 +947,22 @@ int yl,yh; const psImage *image = interp->image; // Image to interpolate - // const psImage *mask = interp->mask; // Mask to interpolate - const psImage *variance = interp->variance; // Variance to interpolate xl = floor(x); if (xl < 0) { xl = 0; } if (xl >= image->numCols) {xl = image->numCols - 1; } + //{ *imageValue = NAN; return PS_INTERPOLATE_STATUS_GOOD; } xh = xl + 1; if (xh >= image->numCols) {xh = image->numCols - 1; } + //{ *imageValue = NAN; return PS_INTERPOLATE_STATUS_GOOD; } + // yl = floor(y); if (yl < 0) { yl = 0; } if (yl >= image->numRows) {yl = image->numRows - 1; } + //{ *imageValue = NAN; return PS_INTERPOLATE_STATUS_GOOD; } + // yh = yl + 1; if (yh >= image->numRows) {yh = image->numRows - 1; } - if (varianceValue && variance) { - u1 = (xh - x) * variance->data.F32[yl][xl] + (x - xl) * variance->data.F32[yl][xh]; - u2 = (xh - x) * variance->data.F32[yh][xl] + (x - xl) * variance->data.F32[yh][xh]; - - *varianceValue = (yh - y) * u1 + (y - yl) * u2; - } -#if 0 - if (maskValue && mask) { - u1 = (xh - x) * mask->data.F32[yl][xl] + (x - xl) * mask->data.F32[yl][xh]; - u2 = (xh - x) * mask->data.F32[yh][xl] + (x - xl) * mask->data.F32[yh][xh]; - - maskValue = (yh - y) * u1 + (y - yl) * u2; - } -#endif + //{ *imageValue = NAN; return PS_INTERPOLATE_STATUS_GOOD; } + // if (imageValue && image) { if (image->data.F32[yl][xl] == image->data.F32[yl][xh]) { Index: /branches/czw_branch/20120906/pswarp/src/pswarp.h =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarp.h (revision 34469) +++ /branches/czw_branch/20120906/pswarp/src/pswarp.h (revision 34470) @@ -94,4 +94,5 @@ pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix); +pswarpMapGrid *pswarpMapGridFromImageForBackground (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); Index: /branches/czw_branch/20120906/pswarp/src/pswarpLoop.c =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarpLoop.c (revision 34469) +++ /branches/czw_branch/20120906/pswarp/src/pswarpLoop.c (revision 34470) @@ -804,5 +804,8 @@ fprintf(stderr,"Transforming Readout!\n"); + + psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", true); pswarpTransformReadout(output, readout, config); + psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", false); if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { Index: /branches/czw_branch/20120906/pswarp/src/pswarpMapGrid.c =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarpMapGrid.c (revision 34469) +++ /branches/czw_branch/20120906/pswarp/src/pswarpMapGrid.c (revision 34470) @@ -52,4 +52,45 @@ return grid; } +/** + * pswarpMapGridFromImageForBackground builds a set (a grid) of locally-linear maps which convert the source + * 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 every source image pixel is guaranteed to have a map. + */ +pswarpMapGrid *pswarpMapGridFromImageForBackground (pmReadout *dest, pmReadout *src, int nXpix, int nYpix) { + + int i, ni; + int j, nj; + + // start counting from the center of the superpixels + double xMin = 0.5*nXpix; + double yMin = 0.5*nYpix; + + // the map is defined for coordinates in the image parent frame. + int Nx = src->image->numCols + src->image->col0; + int Ny = src->image->numRows + src->image->row0; + + // always allocate an extra superpixel to handle spillover + int nXpts = (int)(Nx / nXpix) + 1; + int nYpts = (int)(Ny / nYpix) + 1; + + // create the grid of maps + pswarpMapGrid *grid = pswarpMapGridAlloc (nXpts, nYpts); + + // measure the map for the center of each superpixel + for (ni = 0, i = xMin; ni < nXpts; i += nXpix, ni++) { + for (nj = 0, j = yMin; nj < nYpts; j += nYpix, nj++) { + pswarpMapSetLocalModel (grid->maps[ni][nj], dest, src, i, j); +/* grid->maps[ni][nj]->Xo += (13.0 * 400.0 - 4846.0) / 2.0; */ +/* grid->maps[ni][nj]->Yo += (13.0 * 400.0 - 4868.0) / 2.0; */ + } + } + + grid->nXpix = nXpix; + grid->nYpix = nYpix; + grid->xMin = xMin; + grid->yMin = yMin; + return grid; +} /** Index: /branches/czw_branch/20120906/pswarp/src/pswarpTransformReadout.c =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarpTransformReadout.c (revision 34469) +++ /branches/czw_branch/20120906/pswarp/src/pswarpTransformReadout.c (revision 34470) @@ -60,6 +60,11 @@ // pswarpMapGridFromImage builds a set of locally-linear maps which convert the // output coordinates to input coordinates - pswarpMapGrid *grid = pswarpMapGridFromImage(input, output, nGridX, nGridY); - + pswarpMapGrid *grid; + if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND_WARPING")) { + grid = pswarpMapGridFromImageForBackground(input, output, nGridX, nGridY); + } + else { + grid = pswarpMapGridFromImage(input, output, nGridX, nGridY); + } // XXX optionally modify the grid based on this result and force the maxError < XXX double maxError = pswarpMapGridMaxError(grid); // Maximum (positional) error from using grid Index: /branches/czw_branch/20120906/pswarp/src/pswarpTransformTile.c =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarpTransformTile.c (revision 34469) +++ /branches/czw_branch/20120906/pswarp/src/pswarpTransformTile.c (revision 34470) @@ -86,51 +86,31 @@ for (int y = yMin; y < yMax; y++) { for (int x = xMin; x < xMax; x++) { -/* if ((x == 931)||(x == 717)||(x == 924)) { */ -/* /\* if ((y == 1348)||(y == 1249)||(y == 1081)) { *\/ */ -/* /\* fprintf(stderr,"BEcause: %d %d [%d %d %d %d]\n",x,y,yMin,yMax,xMin,xMax); *\/ */ -/* /\* } *\/ */ -/* } */ - if (((x == 931)&&(y == 1348))|| - ((x == 717)&&(y == 1249))|| - ((x == 924)&&(y == 1081))) { - psTrace("pswarp",3,"TT: A\n"); - psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax); - - } - // Only transform those pixels requested if (region && region->data.PS_TYPE_IMAGE_MASK_DATA[y][x]) { continue; } - if (((x == 931)&&(y == 1348))|| - ((x == 717)&&(y == 1249))|| - ((x == 924)&&(y == 1081))) { - psTrace("pswarp",3,"TT: B\n"); - psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax); - - } - // pswarpMapApply converts the output coordinate (x,y) to the input coordinate. // both are in the parent frames of the input and output images. double xIn, yIn; // Input pixel coordinates pswarpMapApply(&xIn, &yIn, map, x + 0.5, y + 0.5); - if (((x == 931)&&(y == 1348))|| - ((x == 717)&&(y == 1249))|| - ((x == 924)&&(y == 1081))) { - psTrace("pswarp",3,"TT: Ct\n"); - psTrace("pswarp",3,"TT: %d %d %d %d %g %g %d %d\n",xMin,xMax,yMin,yMax,xIn,yIn,inNumCols,inNumRows); - } + if (args->interp->mode == 8) { + xIn += 177.0 / 400.0; + yIn += 166.0 / 400.0; + if ((xIn > inNumCols - 177.0/ 400.0)|| + (yIn > inNumRows - 166.0/ 400.0)|| + (xIn < 177.0/ 400.0)|| + (yIn < 166.0/ 400.0) + + ) { + fprintf(stderr,"\n"); + continue; + } + } + if (xIn < 0 || xIn >= inNumCols || yIn < 0 || yIn >= inNumRows) { continue; } - if (((x == 931)&&(y == 1348))|| - ((x == 717)&&(y == 1249))|| - ((x == 924)&&(y == 1081))) { - psTrace("pswarp",3,"TT: C\n"); - psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax); - - } // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates @@ -144,8 +124,17 @@ int xOut = x - outCol0, yOut = y - outRow0; ///< Position on output image - if ((xIn >= args->interp->image->numCols - 1)|| - (yIn >= args->interp->image->numRows - 1)) { - imageValue = NAN; + if (((xOut % 500) == 0)&&((yOut % 500) == 0)) { + fprintf(stderr,"%d %d %g %g %d %d (%d %d %d %d) %g %g %g\n", + x,y,xIn,yIn,xOut,yOut,xMin,xMax,yMin,yMax,imageValue, + (4846 - 400 * (12)) / 2.0, + (4868 - 400 * (12)) / 2.0 + + ); } + +/* if ((xIn >= args->interp->image->numCols - 1)|| */ +/* (yIn >= args->interp->image->numRows - 1)) { */ +/* imageValue = NAN; */ +/* } */ if (outImageData) { outImageData[yOut][xOut] = imageValue * jacobian; @@ -159,43 +148,4 @@ goodPixels++; - - if ((x == 931)||(x == 717)||(x == 924)) { - if ((y == 1348)||(y == 1249)||(y == 1081)) { - psTrace("pswarp",3,"TT: %d %d %g %g %d %d %g %g (%g) %ld %g %g %g %g\n",x,y,xIn,yIn,xOut,yOut,imageValue,jacobian, - args->input->image->data.F32[(int) yIn][(int) xIn],(long) args->interp, - 0.0,0.0,0.0,0.0); -/* inImage->data.F32[(int) yIn][(int) xIn], */ -/* inImage->data.F32[(int) yIn+1][(int) xIn], */ -/* inImage->data.F32[(int) yIn+1][(int) xIn+1], */ -/* inImage->data.F32[(int) yIn][(int) xIn+1]); */ - psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax); - } - } - - if (((xOut == 931)&&(yOut == 1348))|| - ((xOut == 717)&&(yOut == 1249))|| - ((xOut == 924)&&(yOut == 1081))) { - psTrace("pswarp",3,"TT: %d %d %g %g %d %d %g %g (%g) %ld %g %g %g %g\n",x,y,xIn,yIn,xOut,yOut,imageValue,jacobian, - args->input->image->data.F32[(int) yIn][(int) xIn],(long) args->interp, - 0.0,0.0,0.0,0.0); -/* inImage->data.F32[(int) yIn][(int) xIn], */ -/* inImage->data.F32[(int) yIn+1][(int) xIn], */ -/* inImage->data.F32[(int) yIn+1][(int) xIn+1], */ -/* inImage->data.F32[(int) yIn][(int) xIn+1]); */ - psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax); - - } - if ((xOut == 443)&&(yOut == 659)) { - psTrace("pswarp",3,"TT: %d %d %g %g %d %d %g %g (%g) %ld\n",x,y,xIn,yIn,xOut,yOut,imageValue,jacobian, - args->input->image->data.F32[(int) yIn][(int) xIn],(long) args->interp); - psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax); - - } - if ((xOut == 524)&&(yOut == 576)) { - psTrace("pswarp",3,"TT: %d %d %g %g %d %d %g %g (%g) %ld\n",x,y,xIn,yIn,xOut,yOut,imageValue,jacobian, - args->input->image->data.F32[(int) yIn][(int) xIn],(long) args->interp); - psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax); - - } } }