IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 16, 2008, 4:38:46 PM (17 years ago)
Author:
eugene
Message:

fixed psImageUnbinPixel and psImageUnbin so they are consistent and correct at the edges

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/imageops/psImageUnbin.c

    r14925 r20766  
    77 *  @author Eugene Magnier, IfA
    88 *
    9  *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2007-09-20 23:55:50 $
     9 *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2008-11-17 02:38:43 $
    1111 *
    1212 *  Copyright 2007 Institute for Astronomy, University of Hawaii
    … …  
    2424#include "psImageBinning.h"
    2525#include "psImageUnbin.h"
     26
     27# define UNBIN_TOL 1e-4
    2628
    2729// interpolate from model to background (~bresenham linear interpolation)
    … …  
    6769            // (Xs,Ys) : (Xe,Ye) : binned pixel centers in unbinned coords
    6870            // corresponding to (Ix,Iy), (Ix+1,Iy+1)
    69             int Xs = PS_MAX (0, PS_MIN (Nx, (Ix + 0.5)*DX - dx));
    70             int Ys = PS_MAX (0, PS_MIN (Ny, (Iy + 0.5)*DY - dy));
     71            // XXX should this be "+ dx" and + dy?
     72            int Xs = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, Ix + 0.5)));
     73            int Ys = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, Iy + 0.5)));
    7174            int Xe = PS_MAX (0, PS_MIN (Nx, Xs + DX));
    7275            int Ye = PS_MAX (0, PS_MIN (Ny, Ys + DY));
    7376
    7477            for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
    75                 float Vxs = (V10 - V00)*(iy - Ys) / DY + V00;
    76                 float Vxe = (V11 - V01)*(iy - Ys) / DY + V01;
     78                float dY = (iy - Ys) / (float) DY;
     79                float rY = 1.0 - dY;
     80                float Vxs = V10*dY + V00*rY;
     81                float Vxe = V11*dY + V01*rY;
     82
     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);
     102
    77103                float dV = (Vxe - Vxs) / DX;
    78104                float V  = Vxs;
    79105                for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
    80106                    vOut[iy][ix] = V;
    81                     //assert (fabs(V - psImageUnbinPixel(ix, iy, in, DX, DY, dx, dy)) < 1e-5*fabs(V));
     107                    // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
    82108                    V += dV;
    83109                }
    … …  
    87113
    88114    // side pixels
    89     int Xs = PS_MAX (0, PS_MIN (Nx, ( 0 + 0.5)*DX - dx));
    90     int Xe = PS_MAX (0, PS_MIN (Nx, (nx - 0.5)*DX - dx));
     115    int Xs = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, 0  + 0.5)));
     116    int Xe = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, nx - 0.5)));
    91117    for (int Iy = 0; Iy < ny - 1; Iy++) {
    92118
    93         int Ys = PS_MAX (0, PS_MIN (Ny, (Iy + 0.5)*DY - dy));
     119        int Ys = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, Iy + 0.5)));
    94120        int Ye = PS_MAX (0, PS_MIN (Ny, Ys + DY));
    95121
    … …  
    102128            for (int ix = 0; ix < Xs; ix++) {
    103129                vOut[iy][ix] = V;
     130                // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
    104131            }
    105132            V += dV;
    … …  
    114141            for (int ix = Xe; ix < Nx; ix++) {
    115142                vOut[iy][ix] = V;
     143                // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
    116144            }
    117145            V += dV;
    … …  
    120148
    121149    // top and bottom pixels
    122     int Ys = PS_MAX (0, PS_MIN (Ny, ( 0 + 0.5)*DY - dy));
    123     int Ye = PS_MAX (0, PS_MIN (Ny, (ny - 0.5)*DY - dy));
     150    int Ys = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, 0  + 0.5)));
     151    int Ye = PS_MAX (0, PS_MIN (Ny, psImageBinningGetFineY(binning, ny - 0.5)));
    124152    for (int Ix = 0; Ix < nx - 1; Ix++) {
    125153
    126         int Xs = PS_MAX (0, PS_MIN (Nx, (Ix + 0.5)*DX - dx));
     154        int Xs = PS_MAX (0, PS_MIN (Nx, psImageBinningGetFineX(binning, Ix + 0.5)));
    127155        int Xe = PS_MAX (0, PS_MIN (Nx, Xs + DX));
    128156
    … …  
    135163            for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
    136164                vOut[iy][ix] = V;
     165                // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
    137166                V += dV;
    138167            }
    … …  
    147176            for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
    148177                vOut[iy][ix] = V;
     178                // assert (fabs(V - psImageUnbinPixel(ix, iy, in, binning)) < UNBIN_TOL*fabs(V));
    149179                V += dV;
    150180            }
    … …  
    157187        float V;
    158188        // center of last pixel
    159         int Xs = PS_MAX (0, PS_MIN (Nx, ( 0 + 0.5)*DX - dx));
    160         int Xe = PS_MAX (0, PS_MIN (Nx, (nx - 0.5)*DX - dx));
    161         int Ys = PS_MAX (0, PS_MIN (Ny, ( 0 + 0.5)*DY - dy));
    162         int Ye = PS_MAX (0, PS_MIN (Ny, (ny - 0.5)*DY - dy));
     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)));
    163193
    164194        // 0,0
    165195        V = vIn[0][0];
     196        // assert (fabs(V - psImageUnbinPixel(0, 0, in, binning)) < UNBIN_TOL*fabs(V));
     197
    166198        for (int iy = 0; iy < Ys; iy++)
    167199        {
    … …  
    172204        // Nx,0
    173205        V = vIn[0][nx-1];
     206        // assert (fabs(V - psImageUnbinPixel(Nx-1, 0, in, binning)) < UNBIN_TOL*fabs(V));
     207
    174208        for (int iy = 0; iy < Ys; iy++)
    175209        {
    … …  
    180214        // 0,Ny
    181215        V = vIn[ny-1][0];
     216        // assert (fabs(V - psImageUnbinPixel(0, Ny-1, in, binning)) < UNBIN_TOL*fabs(V));
     217
    182218        for (int iy = Ye; iy < Ny; iy++)
    183219        {
    … …  
    188224        // Nx,Ny
    189225        V = vIn[ny-1][nx-1];
     226        // assert (fabs(V - psImageUnbinPixel(Nx-1, Ny-1, in, binning)) < UNBIN_TOL*fabs(V));
     227
    190228        for (int iy = Ye; iy < Ny; iy++)
    191229        {
    … …  
    202240/*
    203241 * Get the value of a single unbinned pixel from the binned representation
    204  *
     242 */
     243
     244double psImageUnbinPixel(const double xFine, const double yFine, // desired Unbinned point (parent coords)
     245                         const psImage *in, // binned image
     246                         const psImageBinning *binning)   //!< Overhang
     247{
     248    PS_ASSERT_IMAGE_NON_NULL(in, NAN);
     249    assert (in->type.type == PS_TYPE_F32);
     250
     251    const float xRuff = psImageBinningGetRuffX (binning, xFine);
     252    const float yRuff = psImageBinningGetRuffY (binning, yFine);
     253
     254    const double value = psImageInterpolatePixelBilinear (xRuff, yRuff, in);
     255
     256    return value;
     257}
     258
     259// fast & simple API to interpolate to a subpixel position using bilinear interpolation
     260// x,y in parent image coordinates (pixel centers at 0.5, 0.5)
     261double psImageInterpolatePixelBilinear (const double xIn, const double yIn, const psImage *in) {
     262
     263    PS_ASSERT_PTR_NON_NULL(in, PS_ERR_BAD_PARAMETER_VALUE);
     264    assert (in->type.type == PS_TYPE_F32);
     265
     266    const double x = xIn - in->col0;
     267    const double y = yIn - in->row0;
     268
     269    // allow extrapolation one pixel beyong edge of valid pixels, but no further
     270    // (this allows the nXskip,nYskip boundary areas to be used as well)
     271    if ((x < -1) || (x > in->numCols) || (y < -1) || (y > in->numRows)) {
     272        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Point (%f,%f) lies outside binned image", x, y);
     273        return NAN;
     274    }
     275
     276    // limiting cases: Nx == 1 and/or Ny == 1
     277
     278    // if we have a single pixel, there is no spatial information
     279    if ((in->numCols == 1) && (in->numRows == 1)) {
     280        const double value = in->data.F32[0][0];
     281        return value;
     282    }
     283
     284    // handle edge cases with extrapolation
     285
     286    const int ix = x - 0.5; // index of reference pixel
     287    const int iy = y - 0.5; // index of reference pixel
     288
     289    // do numCols,Rows first so we are never < 0
     290    const int Xs = PS_MAX (PS_MIN (ix, in->numCols - 2), 0);
     291    const int Ys = PS_MAX (PS_MIN (iy, in->numRows - 2), 0);
     292
     293    const int Xe = Xs + 1;
     294    const int Ye = Ys + 1;
     295
     296    // dx,dy range from 0.0 to 1.0 for interpolated pixels, and -0.5 to 1.5 for extrapolation
     297    const double dx = x - 0.5 - Xs;
     298    const double dy = y - 0.5 - Ys;
     299
     300    const double rx = 1.0 - dx;
     301    const double ry = 1.0 - dy;
     302
     303    // if Nx == 1, we have no x-dir spatial information
     304    if (in->numCols == 1) {
     305        double V0 = in->data.F32[Ys][Xs];
     306        double V1 = in->data.F32[Ye][Xs];
     307
     308        const double value = V0*ry + V1*dy;
     309        return value;
     310    }   
     311
     312    // if Ny == 1, we have no y-dir spatial information
     313    if (in->numRows == 1) {
     314        double V0 = in->data.F32[Ys][Xs];
     315        double V1 = in->data.F32[Ys][Xe];
     316
     317        const double value = V0*rx + V1*dx;
     318        return value;
     319    }   
     320
     321    // Vxy
     322    double V00 = in->data.F32[Ys][Xs];
     323    double V01 = in->data.F32[Ye][Xs];
     324    double V10 = in->data.F32[Ys][Xe];
     325    double V11 = in->data.F32[Ye][Xe];
     326
     327    double value;
     328
     329    // corners
     330    if ((dx < 0.0) && (dy < 0.0)) {
     331        return V00;
     332    }   
     333    if ((dx > 1.0) && (dy < 0.0)) {
     334        return V10;
     335    }   
     336    if ((dx < 0.0) && (dy > 1.0)) {
     337        return V01;
     338    }   
     339    if ((dx > 1.0) && (dy > 1.0)) {
     340        return V11;
     341    }   
     342
     343    // sides
     344    if (dx < 0.0) {
     345        value = V00*ry + V01*dy;
     346        return value;
     347    }   
     348    if (dy < 0.0) {
     349        value = V00*rx + V10*dx;
     350        return value;
     351    }   
     352    if (dx > 1.0) {
     353        value = V10*ry + V11*dy;
     354        return value;
     355    }   
     356    if (dy > 1.0) {
     357        value = V01*rx + V11*dx;
     358        return value;
     359    }   
     360
     361    // bilinear interpolation
     362    value = V00*rx*ry + V10*dx*ry + V01*rx*dy + V11*dx*dy;
     363    return value;
     364}
     365
     366# if (0)
     367/* Old version of this function
    205368 * N.b. This code only works for the central part of the image; the edge
    206369 * cases should be added
    207 
    208370 * XXXX this is bilinear interpolation, but written sub-optimally
    209371 * XXXX this function should be taking float input coordinates!!!
    … …  
    253415}
    254416
    255 double psImageUnbinPixel_V2(const double xFine, const double yFine, // desired Unbinned point (parent coords)
    256                             const psImage *in, // binned image
    257                             const psImageBinning *binning)   //!< Overhang
    258 {
    259     PS_ASSERT_IMAGE_NON_NULL(in, NAN);
    260     assert (in->type.type == PS_TYPE_F32);
    261 
    262     const float xRuff = psImageBinningGetRuffX (binning, xFine);
    263     const float yRuff = psImageBinningGetRuffY (binning, yFine);
    264 
    265     const double value = psImageInterpolatePixelBilinear (xRuff, yRuff, in);
    266 
    267     return value;
    268 }
    269 
    270 // fast & simple API to interpolate to a subpixel position using bilinear interpolation
    271 // x,y in parent image coordinates (pixel centers at 0.5, 0.5)
    272 double psImageInterpolatePixelBilinear (const double xIn, const double yIn, const psImage *in) {
    273 
    274     PS_ASSERT_PTR_NON_NULL(in, PS_ERR_BAD_PARAMETER_VALUE);
    275     assert (in->type.type == PS_TYPE_F32);
    276 
    277     const double x = xIn - in->col0;
    278     const double y = yIn - in->row0;
    279 
    280     // allow extrapolation to edge of valid pixels, but not beyond
    281     if ((x < 0) || (x >= in->numCols) || (y < 0) || (y >= in->numRows)) {
    282         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Point (%f,%f) lies outside binned image", x, y);
    283         return NAN;
    284     }
    285 
    286     // limiting cases: Nx == 1 and/or Ny == 1
    287 
    288     // if we have a single pixel, there is no spatial information
    289     if ((in->numCols == 1) && (in->numRows == 1)) {
    290         const double value = in->data.F32[0][0];
    291         return value;
    292     }
    293 
    294     // handle edge cases with extrapolation
    295 
    296     const int ix = x - 0.5; // index of reference pixel
    297     const int iy = y - 0.5; // index of reference pixel
    298 
    299     // do numCols,Rows first so we are never < 0
    300     const int Xs = PS_MAX (PS_MIN (ix, in->numCols - 2), 0);
    301     const int Ys = PS_MAX (PS_MIN (iy, in->numRows - 2), 0);
    302 
    303     const int Xe = Xs + 1;
    304     const int Ye = Ys + 1;
    305 
    306     // dx,dy range from 0.0 to 1.0 for interpolated pixels, and -0.5 to 1.5 for extrapolation
    307     const double dx = x - 0.5 - Xs;
    308     const double dy = y - 0.5 - Ys;
    309 
    310     const double rx = 1.0 - dx;
    311     const double ry = 1.0 - dy;
    312 
    313     // if Nx == 1, we have no x-dir spatial information
    314     if (in->numCols == 1) {
    315         double V0 = in->data.F32[Ys][Xs];
    316         double V1 = in->data.F32[Ye][Xs];
    317 
    318         const double value = V0*ry + V1*dy;
    319         return value;
    320     }   
    321 
    322     // if Ny == 1, we have no y-dir spatial information
    323     if (in->numRows == 1) {
    324         double V0 = in->data.F32[Ys][Xs];
    325         double V1 = in->data.F32[Ys][Xe];
    326 
    327         const double value = V0*rx + V1*dx;
    328         return value;
    329     }   
    330 
    331     // Vxy
    332     double V00 = in->data.F32[Ys][Xs];
    333     double V10 = in->data.F32[Ys][Xe];
    334     double V01 = in->data.F32[Ye][Xs];
    335     double V11 = in->data.F32[Ye][Xe];
    336 
    337     // bilinear interpolation
    338     const double value = V00*rx*ry + V10*dx*ry + V01*rx*dy + V11*dx*dy;
    339 
    340     return value;
    341 }
    342 
    343    
     417# endif
     418
Note: See TracChangeset for help on using the changeset viewer.