IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 18, 2009, 4:28:43 PM (17 years ago)
Author:
Paul Price
Message:

Doubles should use lf.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/psLib/src/imageops/psImageUnbin.c

    r21040 r24497  
    3030// DX, DY are the binning factor
    3131// dx, dy is the distance is high-res pixels to the 0,0 corner of the first
    32 // binned pixel. 
     32// binned pixel.
    3333// XXX check that this is still consistent with psphotImageMedian...
    3434psImage *psImageUnbin(psImage *out, const psImage *in, const psImageBinning *binning)
     
    6969            // (Xs,Ys) : (Xe,Ye) : binned pixel centers in unbinned coords
    7070            // corresponding to (Ix,Iy), (Ix+1,Iy+1)
    71             // XXX should this be "+ dx" and + dy?
     71            // XXX should this be "+ dx" and + dy?
    7272            int Xs = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, Ix + 0.5)));
    7373            int Ys = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, Iy + 0.5)));
     
    7676
    7777            for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
    78                 float dY = (iy - Ys) / (float) DY;
    79                 float rY = 1.0 - dY;
     78                float dY = (iy - Ys) / (float) DY;
     79                float rY = 1.0 - dY;
    8080                float Vxs = V10*dY + V00*rY;
    8181                float Vxe = V11*dY + V01*rY;
    8282
    83                 // Vxs = (V10 - V00)*(iy - Ys) / DY + V00;                                                                                                       
    84                 // Vxe = (V11 - V01)*(iy - Ys) / DY + V01;                                                                                                       
    85 
    86                 // dVx = Vxs_1 - Vxs_1
    87                 // dVx = (V10*dY_1 + V00*rY_1) - (V10*dY_0 + V00*rY_0);
    88                 // dY_0 = (iy - Ys)/DY     = iy/DY - Ys/DY;
    89                 // dY_1 = (iy + 1 - Ys)/DY = iy/DY - Ys/DY + 1/DY;
    90                 // rY_0 = 1 - dY_0;
    91                 // rY_1 = 1 - dY_1;
    92                 // dVx = V10*(ddY) + V00*(drY);
    93                 // ddY = 1/DY;
    94                 // drY = -1/DY;
    95 
    96                 // dVxs = (V10 - V00)/DY;
    97                 // dVxe = (V11 - V01)/DY;
    98 
    99                 // ddV = (Vxe_1 - Vxs_1)/DX - (Vxe_0 - Vxs_0)/DX;
    100                 // ddV = (Vxe_1 - Vxe_0)/DX - (Vxs_1 - Vxs_0)/DX;
    101                 // ddV = (V11 - V01 - V10 + V00)/(DX*DY);
     83                // Vxs = (V10 - V00)*(iy - Ys) / DY + V00;
     84                // Vxe = (V11 - V01)*(iy - Ys) / DY + V01;
     85
     86                // dVx = Vxs_1 - Vxs_1
     87                // dVx = (V10*dY_1 + V00*rY_1) - (V10*dY_0 + V00*rY_0);
     88                // dY_0 = (iy - Ys)/DY     = iy/DY - Ys/DY;
     89                // dY_1 = (iy + 1 - Ys)/DY = iy/DY - Ys/DY + 1/DY;
     90                // rY_0 = 1 - dY_0;
     91                // rY_1 = 1 - dY_1;
     92                // dVx = V10*(ddY) + V00*(drY);
     93                // ddY = 1/DY;
     94                // drY = -1/DY;
     95
     96                // dVxs = (V10 - V00)/DY;
     97                // dVxe = (V11 - V01)/DY;
     98
     99                // ddV = (Vxe_1 - Vxs_1)/DX - (Vxe_0 - Vxs_0)/DX;
     100                // ddV = (Vxe_1 - Vxe_0)/DX - (Vxs_1 - Vxs_0)/DX;
     101                // ddV = (V11 - V01 - V10 + V00)/(DX*DY);
    102102
    103103                float dV = (Vxe - Vxs) / DX;
     
    128128            for (int ix = 0; ix < Xs; ix++) {
    129129                vOut[iy][ix] = V;
    130                 // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
     130                // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
    131131            }
    132132            V += dV;
     
    141141            for (int ix = Xe; ix < Nx; ix++) {
    142142                vOut[iy][ix] = V;
    143                 // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
     143                // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
    144144            }
    145145            V += dV;
     
    163163            for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
    164164                vOut[iy][ix] = V;
    165                 // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
     165                // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
    166166                V += dV;
    167167            }
     
    176176            for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
    177177                vOut[iy][ix] = V;
    178                 // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
     178                // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
    179179                V += dV;
    180180            }
     
    186186    {
    187187        float V;
    188         // center of last pixel
    189         int Xs = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, 0  + 0.5)));
    190         int Xe = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, nx - 0.5)));
    191         int Ys = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, 0  + 0.5)));
    192         int Ye = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, ny - 0.5)));
     188        // center of last pixel
     189        int Xs = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, 0  + 0.5)));
     190        int Xe = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, nx - 0.5)));
     191        int Ys = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, 0  + 0.5)));
     192        int Ye = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, ny - 0.5)));
    193193
    194194        // 0,0
    195195        V = vIn[0][0];
    196         // assert (fabs(V - psImageUnbinPixel(0, 0, in, binning)) < UNBIN_TOL*fabs(V));
     196        // assert (fabs(V - psImageUnbinPixel(0, 0, in, binning)) < UNBIN_TOL*fabs(V));
    197197
    198198        for (int iy = 0; iy < Ys; iy++)
     
    204204        // Nx,0
    205205        V = vIn[0][nx-1];
    206         // assert (fabs(V - psImageUnbinPixel(Nx-1, 0, in, binning)) < UNBIN_TOL*fabs(V));
     206        // assert (fabs(V - psImageUnbinPixel(Nx-1, 0, in, binning)) < UNBIN_TOL*fabs(V));
    207207
    208208        for (int iy = 0; iy < Ys; iy++)
     
    214214        // 0,Ny
    215215        V = vIn[ny-1][0];
    216         // assert (fabs(V - psImageUnbinPixel(0, Ny-1, in, binning)) < UNBIN_TOL*fabs(V));
     216        // assert (fabs(V - psImageUnbinPixel(0, Ny-1, in, binning)) < UNBIN_TOL*fabs(V));
    217217
    218218        for (int iy = Ye; iy < Ny; iy++)
     
    224224        // Nx,Ny
    225225        V = vIn[ny-1][nx-1];
    226         // assert (fabs(V - psImageUnbinPixel(Nx-1, Ny-1, in, binning)) < UNBIN_TOL*fabs(V));
     226        // assert (fabs(V - psImageUnbinPixel(Nx-1, Ny-1, in, binning)) < UNBIN_TOL*fabs(V));
    227227
    228228        for (int iy = Ye; iy < Ny; iy++)
     
    243243
    244244double psImageUnbinPixel(const double xFine, const double yFine, // desired Unbinned point (parent coords)
    245                         const psImage *in, // binned image
    246                         const psImageBinning *binning)   //!< Overhang
     245                        const psImage *in, // binned image
     246                        const psImageBinning *binning)   //!< Overhang
    247247{
    248248    PS_ASSERT_IMAGE_NON_NULL(in, NAN);
     
    273273
    274274    if ((x < -nXedge) || (x > in->numCols + nXedge) || (y < -nYedge) || (y > in->numRows + nYedge)) {
    275         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Point (%f,%f) lies outside binned image", x, y);
     275        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Point (%lf,%lf) lies outside binned image", x, y);
    276276        return NAN;
    277277    }
     
    281281    // if we have a single pixel, there is no spatial information
    282282    if ((in->numCols == 1) && (in->numRows == 1)) {
    283         const double value = in->data.F32[0][0];
    284         return value;
     283        const double value = in->data.F32[0][0];
     284        return value;
    285285    }
    286286
     
    306306    // if Nx == 1, we have no x-dir spatial information
    307307    if (in->numCols == 1) {
    308         double V0 = in->data.F32[Ys][Xs];
    309         double V1 = in->data.F32[Ye][Xs];
    310 
    311         const double value = V0*ry + V1*dy;
    312         return value;
    313     }   
     308        double V0 = in->data.F32[Ys][Xs];
     309        double V1 = in->data.F32[Ye][Xs];
     310
     311        const double value = V0*ry + V1*dy;
     312        return value;
     313    }
    314314
    315315    // if Ny == 1, we have no y-dir spatial information
    316316    if (in->numRows == 1) {
    317         double V0 = in->data.F32[Ys][Xs];
    318         double V1 = in->data.F32[Ys][Xe];
    319 
    320         const double value = V0*rx + V1*dx;
    321         return value;
    322     }   
    323 
    324     // Vxy 
     317        double V0 = in->data.F32[Ys][Xs];
     318        double V1 = in->data.F32[Ys][Xe];
     319
     320        const double value = V0*rx + V1*dx;
     321        return value;
     322    }
     323
     324    // Vxy
    325325    double V00 = in->data.F32[Ys][Xs];
    326326    double V01 = in->data.F32[Ye][Xs];
     
    332332    // corners
    333333    if ((dx < 0.0) && (dy < 0.0)) {
    334         return V00;
    335     }   
     334        return V00;
     335    }
    336336    if ((dx > 1.0) && (dy < 0.0)) {
    337         return V10;
    338     }   
     337        return V10;
     338    }
    339339    if ((dx < 0.0) && (dy > 1.0)) {
    340         return V01;
    341     }   
     340        return V01;
     341    }
    342342    if ((dx > 1.0) && (dy > 1.0)) {
    343         return V11;
    344     }   
     343        return V11;
     344    }
    345345
    346346    // sides
    347347    if (dx < 0.0) {
    348         value = V00*ry + V01*dy;
    349         return value;
    350     }   
     348        value = V00*ry + V01*dy;
     349        return value;
     350    }
    351351    if (dy < 0.0) {
    352         value = V00*rx + V10*dx;
    353         return value;
    354     }   
     352        value = V00*rx + V10*dx;
     353        return value;
     354    }
    355355    if (dx > 1.0) {
    356         value = V10*ry + V11*dy;
    357         return value;
    358     }   
     356        value = V10*ry + V11*dy;
     357        return value;
     358    }
    359359    if (dy > 1.0) {
    360         value = V01*rx + V11*dx;
    361         return value;
    362     }   
     360        value = V01*rx + V11*dx;
     361        return value;
     362    }
    363363
    364364    // bilinear interpolation
Note: See TracChangeset for help on using the changeset viewer.