Index: trunk/psLib/src/imageops/psImageUnbin.c =================================================================== --- trunk/psLib/src/imageops/psImageUnbin.c (revision 14925) +++ trunk/psLib/src/imageops/psImageUnbin.c (revision 20766) @@ -7,6 +7,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-09-20 23:55:50 $ + * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-11-17 02:38:43 $ * * Copyright 2007 Institute for Astronomy, University of Hawaii @@ -24,4 +24,6 @@ #include "psImageBinning.h" #include "psImageUnbin.h" + +# define UNBIN_TOL 1e-4 // interpolate from model to background (~bresenham linear interpolation) @@ -67,17 +69,41 @@ // (Xs,Ys) : (Xe,Ye) : binned pixel centers in unbinned coords // corresponding to (Ix,Iy), (Ix+1,Iy+1) - int Xs = PS_MAX (0, PS_MIN (Nx, (Ix + 0.5)*DX - dx)); - int Ys = PS_MAX (0, PS_MIN (Ny, (Iy + 0.5)*DY - dy)); + // XXX should this be "+ dx" and + dy? + int Xs = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, Ix + 0.5))); + int Ys = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, Iy + 0.5))); int Xe = PS_MAX (0, PS_MIN (Nx, Xs + DX)); int Ye = PS_MAX (0, PS_MIN (Ny, Ys + DY)); for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) { - float Vxs = (V10 - V00)*(iy - Ys) / DY + V00; - float Vxe = (V11 - V01)*(iy - Ys) / DY + V01; + float dY = (iy - Ys) / (float) DY; + float rY = 1.0 - dY; + float Vxs = V10*dY + V00*rY; + float Vxe = V11*dY + V01*rY; + + // Vxs = (V10 - V00)*(iy - Ys) / DY + V00; + // Vxe = (V11 - V01)*(iy - Ys) / DY + V01; + + // dVx = Vxs_1 - Vxs_1 + // dVx = (V10*dY_1 + V00*rY_1) - (V10*dY_0 + V00*rY_0); + // dY_0 = (iy - Ys)/DY = iy/DY - Ys/DY; + // dY_1 = (iy + 1 - Ys)/DY = iy/DY - Ys/DY + 1/DY; + // rY_0 = 1 - dY_0; + // rY_1 = 1 - dY_1; + // dVx = V10*(ddY) + V00*(drY); + // ddY = 1/DY; + // drY = -1/DY; + + // dVxs = (V10 - V00)/DY; + // dVxe = (V11 - V01)/DY; + + // ddV = (Vxe_1 - Vxs_1)/DX - (Vxe_0 - Vxs_0)/DX; + // ddV = (Vxe_1 - Vxe_0)/DX - (Vxs_1 - Vxs_0)/DX; + // ddV = (V11 - V01 - V10 + V00)/(DX*DY); + float dV = (Vxe - Vxs) / DX; float V = Vxs; for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) { vOut[iy][ix] = V; - //assert (fabs(V - psImageUnbinPixel(ix, iy, in, DX, DY, dx, dy)) < 1e-5*fabs(V)); + // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V)); V += dV; } @@ -87,9 +113,9 @@ // side pixels - int Xs = PS_MAX (0, PS_MIN (Nx, ( 0 + 0.5)*DX - dx)); - int Xe = PS_MAX (0, PS_MIN (Nx, (nx - 0.5)*DX - dx)); + int Xs = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, 0 + 0.5))); + int Xe = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, nx - 0.5))); for (int Iy = 0; Iy < ny - 1; Iy++) { - int Ys = PS_MAX (0, PS_MIN (Ny, (Iy + 0.5)*DY - dy)); + int Ys = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, Iy + 0.5))); int Ye = PS_MAX (0, PS_MIN (Ny, Ys + DY)); @@ -102,4 +128,5 @@ for (int ix = 0; ix < Xs; ix++) { vOut[iy][ix] = V; + // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V)); } V += dV; @@ -114,4 +141,5 @@ for (int ix = Xe; ix < Nx; ix++) { vOut[iy][ix] = V; + // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V)); } V += dV; @@ -120,9 +148,9 @@ // top and bottom pixels - int Ys = PS_MAX (0, PS_MIN (Ny, ( 0 + 0.5)*DY - dy)); - int Ye = PS_MAX (0, PS_MIN (Ny, (ny - 0.5)*DY - dy)); + int Ys = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, 0 + 0.5))); + int Ye = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, ny - 0.5))); for (int Ix = 0; Ix < nx - 1; Ix++) { - int Xs = PS_MAX (0, PS_MIN (Nx, (Ix + 0.5)*DX - dx)); + int Xs = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, Ix + 0.5))); int Xe = PS_MAX (0, PS_MIN (Nx, Xs + DX)); @@ -135,4 +163,5 @@ for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) { vOut[iy][ix] = V; + // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V)); V += dV; } @@ -147,4 +176,5 @@ for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) { vOut[iy][ix] = V; + // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V)); V += dV; } @@ -157,11 +187,13 @@ float V; // center of last pixel - int Xs = PS_MAX (0, PS_MIN (Nx, ( 0 + 0.5)*DX - dx)); - int Xe = PS_MAX (0, PS_MIN (Nx, (nx - 0.5)*DX - dx)); - int Ys = PS_MAX (0, PS_MIN (Ny, ( 0 + 0.5)*DY - dy)); - int Ye = PS_MAX (0, PS_MIN (Ny, (ny - 0.5)*DY - dy)); + int Xs = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, 0 + 0.5))); + int Xe = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, nx - 0.5))); + int Ys = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, 0 + 0.5))); + int Ye = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, ny - 0.5))); // 0,0 V = vIn[0][0]; + // assert (fabs(V - psImageUnbinPixel(0, 0, in, binning)) < UNBIN_TOL*fabs(V)); + for (int iy = 0; iy < Ys; iy++) { @@ -172,4 +204,6 @@ // Nx,0 V = vIn[0][nx-1]; + // assert (fabs(V - psImageUnbinPixel(Nx-1, 0, in, binning)) < UNBIN_TOL*fabs(V)); + for (int iy = 0; iy < Ys; iy++) { @@ -180,4 +214,6 @@ // 0,Ny V = vIn[ny-1][0]; + // assert (fabs(V - psImageUnbinPixel(0, Ny-1, in, binning)) < UNBIN_TOL*fabs(V)); + for (int iy = Ye; iy < Ny; iy++) { @@ -188,4 +224,6 @@ // Nx,Ny V = vIn[ny-1][nx-1]; + // assert (fabs(V - psImageUnbinPixel(Nx-1, Ny-1, in, binning)) < UNBIN_TOL*fabs(V)); + for (int iy = Ye; iy < Ny; iy++) { @@ -202,8 +240,132 @@ /* * Get the value of a single unbinned pixel from the binned representation - * + */ + +double psImageUnbinPixel(const double xFine, const double yFine, // desired Unbinned point (parent coords) + const psImage *in, // binned image + const psImageBinning *binning) //!< Overhang +{ + PS_ASSERT_IMAGE_NON_NULL(in, NAN); + assert (in->type.type == PS_TYPE_F32); + + const float xRuff = psImageBinningGetRuffX (binning, xFine); + const float yRuff = psImageBinningGetRuffY (binning, yFine); + + const double value = psImageInterpolatePixelBilinear (xRuff, yRuff, in); + + return value; +} + +// fast & simple API to interpolate to a subpixel position using bilinear interpolation +// x,y in parent image coordinates (pixel centers at 0.5, 0.5) +double psImageInterpolatePixelBilinear (const double xIn, const double yIn, const psImage *in) { + + PS_ASSERT_PTR_NON_NULL(in, PS_ERR_BAD_PARAMETER_VALUE); + assert (in->type.type == PS_TYPE_F32); + + const double x = xIn - in->col0; + const double y = yIn - in->row0; + + // allow extrapolation one pixel beyong edge of valid pixels, but no further + // (this allows the nXskip,nYskip boundary areas to be used as well) + if ((x < -1) || (x > in->numCols) || (y < -1) || (y > in->numRows)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Point (%f,%f) lies outside binned image", x, y); + return NAN; + } + + // limiting cases: Nx == 1 and/or Ny == 1 + + // if we have a single pixel, there is no spatial information + if ((in->numCols == 1) && (in->numRows == 1)) { + const double value = in->data.F32[0][0]; + return value; + } + + // handle edge cases with extrapolation + + const int ix = x - 0.5; // index of reference pixel + const int iy = y - 0.5; // index of reference pixel + + // do numCols,Rows first so we are never < 0 + const int Xs = PS_MAX (PS_MIN (ix, in->numCols - 2), 0); + const int Ys = PS_MAX (PS_MIN (iy, in->numRows - 2), 0); + + const int Xe = Xs + 1; + const int Ye = Ys + 1; + + // dx,dy range from 0.0 to 1.0 for interpolated pixels, and -0.5 to 1.5 for extrapolation + const double dx = x - 0.5 - Xs; + const double dy = y - 0.5 - Ys; + + const double rx = 1.0 - dx; + const double ry = 1.0 - dy; + + // if Nx == 1, we have no x-dir spatial information + if (in->numCols == 1) { + double V0 = in->data.F32[Ys][Xs]; + double V1 = in->data.F32[Ye][Xs]; + + const double value = V0*ry + V1*dy; + return value; + } + + // if Ny == 1, we have no y-dir spatial information + if (in->numRows == 1) { + double V0 = in->data.F32[Ys][Xs]; + double V1 = in->data.F32[Ys][Xe]; + + const double value = V0*rx + V1*dx; + return value; + } + + // Vxy + double V00 = in->data.F32[Ys][Xs]; + double V01 = in->data.F32[Ye][Xs]; + double V10 = in->data.F32[Ys][Xe]; + double V11 = in->data.F32[Ye][Xe]; + + double value; + + // corners + if ((dx < 0.0) && (dy < 0.0)) { + return V00; + } + if ((dx > 1.0) && (dy < 0.0)) { + return V10; + } + if ((dx < 0.0) && (dy > 1.0)) { + return V01; + } + if ((dx > 1.0) && (dy > 1.0)) { + return V11; + } + + // sides + if (dx < 0.0) { + value = V00*ry + V01*dy; + return value; + } + if (dy < 0.0) { + value = V00*rx + V10*dx; + return value; + } + if (dx > 1.0) { + value = V10*ry + V11*dy; + return value; + } + if (dy > 1.0) { + value = V01*rx + V11*dx; + return value; + } + + // bilinear interpolation + value = V00*rx*ry + V10*dx*ry + V01*rx*dy + V11*dx*dy; + return value; +} + +# if (0) +/* Old version of this function * N.b. This code only works for the central part of the image; the edge * cases should be added - * XXXX this is bilinear interpolation, but written sub-optimally * XXXX this function should be taking float input coordinates!!! @@ -253,91 +415,4 @@ } -double psImageUnbinPixel_V2(const double xFine, const double yFine, // desired Unbinned point (parent coords) - const psImage *in, // binned image - const psImageBinning *binning) //!< Overhang -{ - PS_ASSERT_IMAGE_NON_NULL(in, NAN); - assert (in->type.type == PS_TYPE_F32); - - const float xRuff = psImageBinningGetRuffX (binning, xFine); - const float yRuff = psImageBinningGetRuffY (binning, yFine); - - const double value = psImageInterpolatePixelBilinear (xRuff, yRuff, in); - - return value; -} - -// fast & simple API to interpolate to a subpixel position using bilinear interpolation -// x,y in parent image coordinates (pixel centers at 0.5, 0.5) -double psImageInterpolatePixelBilinear (const double xIn, const double yIn, const psImage *in) { - - PS_ASSERT_PTR_NON_NULL(in, PS_ERR_BAD_PARAMETER_VALUE); - assert (in->type.type == PS_TYPE_F32); - - const double x = xIn - in->col0; - const double y = yIn - in->row0; - - // allow extrapolation to edge of valid pixels, but not beyond - if ((x < 0) || (x >= in->numCols) || (y < 0) || (y >= in->numRows)) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Point (%f,%f) lies outside binned image", x, y); - return NAN; - } - - // limiting cases: Nx == 1 and/or Ny == 1 - - // if we have a single pixel, there is no spatial information - if ((in->numCols == 1) && (in->numRows == 1)) { - const double value = in->data.F32[0][0]; - return value; - } - - // handle edge cases with extrapolation - - const int ix = x - 0.5; // index of reference pixel - const int iy = y - 0.5; // index of reference pixel - - // do numCols,Rows first so we are never < 0 - const int Xs = PS_MAX (PS_MIN (ix, in->numCols - 2), 0); - const int Ys = PS_MAX (PS_MIN (iy, in->numRows - 2), 0); - - const int Xe = Xs + 1; - const int Ye = Ys + 1; - - // dx,dy range from 0.0 to 1.0 for interpolated pixels, and -0.5 to 1.5 for extrapolation - const double dx = x - 0.5 - Xs; - const double dy = y - 0.5 - Ys; - - const double rx = 1.0 - dx; - const double ry = 1.0 - dy; - - // if Nx == 1, we have no x-dir spatial information - if (in->numCols == 1) { - double V0 = in->data.F32[Ys][Xs]; - double V1 = in->data.F32[Ye][Xs]; - - const double value = V0*ry + V1*dy; - return value; - } - - // if Ny == 1, we have no y-dir spatial information - if (in->numRows == 1) { - double V0 = in->data.F32[Ys][Xs]; - double V1 = in->data.F32[Ys][Xe]; - - const double value = V0*rx + V1*dx; - return value; - } - - // Vxy - double V00 = in->data.F32[Ys][Xs]; - double V10 = in->data.F32[Ys][Xe]; - double V01 = in->data.F32[Ye][Xs]; - double V11 = in->data.F32[Ye][Xe]; - - // bilinear interpolation - const double value = V00*rx*ry + V10*dx*ry + V01*rx*dy + V11*dx*dy; - - return value; -} - - +# endif +