Index: /trunk/psLib/src/imageops/psImageUnbin.c =================================================================== --- /trunk/psLib/src/imageops/psImageUnbin.c (revision 24496) +++ /trunk/psLib/src/imageops/psImageUnbin.c (revision 24497) @@ -30,5 +30,5 @@ // DX, DY are the binning factor // dx, dy is the distance is high-res pixels to the 0,0 corner of the first -// binned pixel. +// binned pixel. // XXX check that this is still consistent with psphotImageMedian... psImage *psImageUnbin(psImage *out, const psImage *in, const psImageBinning *binning) @@ -69,5 +69,5 @@ // (Xs,Ys) : (Xe,Ye) : binned pixel centers in unbinned coords // corresponding to (Ix,Iy), (Ix+1,Iy+1) - // XXX should this be "+ dx" and + 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))); @@ -76,28 +76,28 @@ for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) { - float dY = (iy - Ys) / (float) DY; - float rY = 1.0 - dY; + 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); + // 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; @@ -128,5 +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)); + // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V)); } V += dV; @@ -141,5 +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)); + // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V)); } V += dV; @@ -163,5 +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)); + // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V)); V += dV; } @@ -176,5 +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)); + // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V)); V += dV; } @@ -186,13 +186,13 @@ { float V; - // center of last pixel - 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))); + // center of last pixel + 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)); + // assert (fabs(V - psImageUnbinPixel(0, 0, in, binning)) < UNBIN_TOL*fabs(V)); for (int iy = 0; iy < Ys; iy++) @@ -204,5 +204,5 @@ // Nx,0 V = vIn[0][nx-1]; - // assert (fabs(V - psImageUnbinPixel(Nx-1, 0, in, binning)) < UNBIN_TOL*fabs(V)); + // assert (fabs(V - psImageUnbinPixel(Nx-1, 0, in, binning)) < UNBIN_TOL*fabs(V)); for (int iy = 0; iy < Ys; iy++) @@ -214,5 +214,5 @@ // 0,Ny V = vIn[ny-1][0]; - // assert (fabs(V - psImageUnbinPixel(0, Ny-1, in, binning)) < UNBIN_TOL*fabs(V)); + // assert (fabs(V - psImageUnbinPixel(0, Ny-1, in, binning)) < UNBIN_TOL*fabs(V)); for (int iy = Ye; iy < Ny; iy++) @@ -224,5 +224,5 @@ // Nx,Ny V = vIn[ny-1][nx-1]; - // assert (fabs(V - psImageUnbinPixel(Nx-1, Ny-1, in, binning)) < UNBIN_TOL*fabs(V)); + // assert (fabs(V - psImageUnbinPixel(Nx-1, Ny-1, in, binning)) < UNBIN_TOL*fabs(V)); for (int iy = Ye; iy < Ny; iy++) @@ -243,6 +243,6 @@ double psImageUnbinPixel(const double xFine, const double yFine, // desired Unbinned point (parent coords) - const psImage *in, // binned image - const psImageBinning *binning) //!< Overhang + const psImage *in, // binned image + const psImageBinning *binning) //!< Overhang { PS_ASSERT_IMAGE_NON_NULL(in, NAN); @@ -273,5 +273,5 @@ if ((x < -nXedge) || (x > in->numCols + nXedge) || (y < -nYedge) || (y > in->numRows + nYedge)) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Point (%f,%f) lies outside binned image", x, y); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Point (%lf,%lf) lies outside binned image", x, y); return NAN; } @@ -281,6 +281,6 @@ // 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; + const double value = in->data.F32[0][0]; + return value; } @@ -306,21 +306,21 @@ // 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; - } + 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 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]; @@ -332,33 +332,33 @@ // corners if ((dx < 0.0) && (dy < 0.0)) { - return V00; - } + return V00; + } if ((dx > 1.0) && (dy < 0.0)) { - return V10; - } + return V10; + } if ((dx < 0.0) && (dy > 1.0)) { - return V01; - } + return V01; + } if ((dx > 1.0) && (dy > 1.0)) { - return V11; - } + return V11; + } // sides if (dx < 0.0) { - value = V00*ry + V01*dy; - return value; - } + value = V00*ry + V01*dy; + return value; + } if (dy < 0.0) { - value = V00*rx + V10*dx; - return value; - } + value = V00*rx + V10*dx; + return value; + } if (dx > 1.0) { - value = V10*ry + V11*dy; - return value; - } + value = V10*ry + V11*dy; + return value; + } if (dy > 1.0) { - value = V01*rx + V11*dx; - return value; - } + value = V01*rx + V11*dx; + return value; + } // bilinear interpolation