IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 12, 2006, 8:59:15 PM (20 years ago)
Author:
eugene
Message:

tests on non-poisson errors, updates to image background, tests of systematic biases, updates to use pmFPAfile paradigm

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotImageMedian.c

    r6753 r6851  
    11# include "psphot.h"
    2 double psImageClippedStats (psImage *image, psImage *mask, psU8 maskValue, double fmin, double fmax);
    3 void fsort (float *value, int N);
    4 
    5 // random number seed to select a fraction of the image pixels
    6 static psRandom *rnd;
    7 
    8 // use no more than MAX_SAMPLE_PIXELS pixels for each median box
    9 static int MAX_SAMPLE_PIXELS;
    102
    113// generate the median in NxN boxes, clipping heavily
    … …  
    157    bool status;
    168    psRegion region;
    17     psImage *model = NULL;
     9    int MAX_SAMPLE_PIXELS;
    1810
    1911    psTimerStart ("psphot");
    2012
     13    // select the appropriate recipe information
     14    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
     15
     16    MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
     17    if (!status) MAX_SAMPLE_PIXELS = 1000;
     18    psImageClippedStatsInit(MAX_SAMPLE_PIXELS);
     19
     20    // subtract this amount extra from the sky
     21    float SKY_BIAS = psMetadataLookupF32 (&status, recipe, "SKY_BIAS");
     22    if (!status) SKY_BIAS = 0;
     23
    2124    // find the currently selected readout
    22     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
    23     pmFPAfile  *input   = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
    24     pmReadout  *readout = pmFPAviewThisReadout (view, input->fpa);
     25    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    2526
    2627    psImage *image = readout->image;
    2728    psImage *mask  = readout->mask;
    28 
    29     rnd = psRandomAlloc (PS_RANDOM_TAUS, 0);
    30     MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
    31     if (!status) MAX_SAMPLE_PIXELS = 1000;
    3229
    3330    // dimensions of input & output image
    … …  
    4946    int ny = (Ny % DY) ? (int)(Ny / DY) + 1 : Ny / DY;
    5047
    51     // select model pixels, from output background model file, or create
    52     model = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKMDL", nx, ny, PS_TYPE_F32);
     48    // select model pixels (from output background model file, or create internal file)
     49    pmReadout *model = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKMDL");
     50    if (model == NULL) {
     51        // select model pixels, from output background model file, or create
     52        model = pmFPAfileCreateInternal (config->files, "PSPHOT.BACKMDL", nx, ny, PS_TYPE_F32);
     53    } else {
     54        // replace the supplied image data with an image of the desired size
     55        psImageRecycle (model->image, nx, ny, PS_TYPE_F32);
     56    }
     57    psF32 **modelData = model->image->data.F32;
    5358
    5459    // measure clipped median for subimages
    … …  
    6368            psImage *submask = psImageSubset (mask, region);
    6469
    65             model->data.F32[iy][ix] = psImageClippedStats (subset, submask, 0xff, 0.25, 0.50);
     70            // XXX the value of the upper and lower cuts probably should be studied...
     71            modelData[iy][ix] = psImageClippedStats (subset, submask, 0xff, 0.25, 0.75) + SKY_BIAS;
    6672
    6773            psFree (subset);
    … …  
    7682    psLogMsg ("psphot", 3, "build median image: %f sec\n", psTimerMark ("psphot"));
    7783
    78     // linear interpolation to full-scale
     84    // XXX temporarily until we get the pmFPAfile output finished
     85    psphotSaveImage (NULL, model->image, "model.fits");
    7986
    8087    // select background pixels, from output background file, or create
    81     psImage *background = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKGND", Nx, Ny, PS_TYPE_F32);
    82 
    83     // XXX this code skips the initial pixels
    84     for (int Iy = 0; Iy < ny-1; Iy ++) {
    85         for (int Ix = 0; Ix < nx-1; Ix ++) {
    86 
    87             float V00 = model->data.F32[Iy+0][Ix+0];
    88             float V01 = model->data.F32[Iy+0][Ix+1];
    89             float V10 = model->data.F32[Iy+1][Ix+0];
    90             float V11 = model->data.F32[Iy+1][Ix+1];
    91 
    92             // a single binned pixel quad
    93             // (Xs,Ys) : (Xe,Ye) : binned pixel centers in unbinned coords
    94             // corresponding to (Ix,Iy), (Ix+1,Iy+1)
    95             int Xs = (Ix + 1)*DX - dx;
    96             int Ys = (Iy + 1)*DY - dy;
    97             int Xe = Xs + DX;
    98             int Ye = Ys + DY;
    99 
    100             for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
    101                 float Vxs = (V10 - V00)*(iy - Ys) / DY + V00;
    102                 float Vxe = (V11 - V01)*(iy - Ys) / DY + V01;
    103                 float dV = (Vxe - Vxs) / DX;
    104                 float V  = Vxs;
    105                 for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
    106                     background->data.F32[iy][ix] = V;
    107                     V += dV;
    108                 }
    109             }
    110         }
     88    pmReadout *background = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKGND");
     89    if (background == NULL) {
     90        background = pmFPAfileCreateInternal (config->files, "PSPHOT.BACKGND", Nx, Ny, PS_TYPE_F32);
     91    } else {
     92        // replace the supplied image data with an image of the desired size
     93        // XXX should not have to do this: it should be allocated correctly
     94        // psImageRecycle (background->image, Nx, Ny, PS_TYPE_F32);
    11195    }
    11296
    113     // side pixels
    114     int Xs = DX - dx;
    115     int Xe = nx*DX - dx;
    116     for (int Iy = 0; Iy < ny - 1; Iy++) {
     97    psF32 **backData = background->image->data.F32;
    11798
    118         int Ys = (Iy + 1)*DY - dy;
    119         int Ye = Ys + DY;
    120 
    121         // leading edge
    122         float V0 = model->data.F32[Iy+0][0];
    123         float V1 = model->data.F32[Iy+1][0];
    124         float dV = (V1 - V0) / DY;
    125         float V = V0;
    126         for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
    127             for (int ix = 0; ix < Xs; ix++) {
    128                 background->data.F32[iy][ix] = V;
    129             }
    130             V += dV;
    131         }
    132 
    133         // trailing edge
    134         V0 = model->data.F32[Iy+0][nx-1];
    135         V1 = model->data.F32[Iy+1][nx-1];
    136         dV = (V1 - V0) / DY;
    137         V = V0;
    138         for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
    139             for (int ix = Xe; ix < Nx; ix++) {
    140                 background->data.F32[iy][ix] = V;
    141             }
    142             V += dV;
    143         }
    144     }
    145 
    146     // top and bottom pixels
    147     int Ys = DY - dy;
    148     int Ye = ny*DY - dy;
    149     for (int Ix = 0; Ix < nx - 1; Ix++) {
    150 
    151         int Xs = (Ix + 1)*DX - dx;
    152         int Xe = Xs + DX;
    153 
    154         // top edge
    155         float V0 = model->data.F32[0][Ix+0];
    156         float V1 = model->data.F32[0][Ix+1];
    157         float dV = (V1 - V0) / DX;
    158         for (int iy = 0; iy < Ys; iy++) {
    159             float V = V0;
    160             for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
    161                 background->data.F32[iy][ix] = V;
    162                 V += dV;
    163             }
    164         }
    165 
    166         // bottom edge
    167         V0 = model->data.F32[ny-1][Ix+0];
    168         V1 = model->data.F32[ny-1][Ix+1];
    169         dV = (V1 - V0) / DX;
    170         for (int iy = Ye; iy < Ny; iy++) {
    171             float V = V0;
    172             for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
    173                 background->data.F32[iy][ix] = V;
    174                 V += dV;
    175             }
    176         }
    177     }
    178 
    179     // the four corners
    180     {
    181         float V;
    182         // 0,0
    183         V = model->data.F32[0][0];
    184         for (int iy = 0; iy < DY - dy; iy++) {
    185             for (int ix = 0; ix < DX - dx; ix++) {
    186                 background->data.F32[iy][ix] = V;
    187             }
    188         }
    189         // Nx,0
    190         V = model->data.F32[0][nx-1];
    191         for (int iy = 0; iy < DY - dy; iy++) {
    192             for (int ix = nx*DX - dx; ix < Nx; ix++) {
    193                 background->data.F32[iy][ix] = V;
    194             }
    195         }
    196         // 0,Ny
    197         V = model->data.F32[ny-1][0];
    198         for (int iy = ny*DY - dy; iy < Ny; iy++) {
    199             for (int ix = 0; ix < DX - dx; ix++) {
    200                 background->data.F32[iy][ix] = V;
    201             }
    202         }
    203         // Nx,Ny
    204         V = model->data.F32[nx-1][ny-1];
    205         for (int iy = ny*DY - dy; iy < Ny; iy++) {
    206             for (int ix = nx*DX - dx; ix < Nx; ix++) {
    207                 background->data.F32[iy][ix] = V;
    208             }
    209         }
    210     }   
     99    // linear interpolation to full-scale
     100    psImageUnbin (background->image, model->image, DX, DY, dx, dy);
    211101    psLogMsg ("psphot", 3, "build resampled image: %f sec\n", psTimerMark ("psphot"));
    212102
    213     psphotSaveImage (NULL, background, "back.fits");
     103    // XXX temporarily until we get the pmFPAfile output finished...
     104    psphotSaveImage (NULL, background->image, "back.fits");
    214105
    215106    // back-sub image pixels, from output background file (don't create if not requested)
    216     psImage *backSub = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKSUB", 0, 0, PS_TYPE_F32);
     107    pmReadout *backSub = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKSUB");
    217108
    218     // subtract the background model
     109    // subtract the background model (save in backSub, if requested)
    219110    for (int j = 0; j < image->numRows; j++) {
    220111        for (int i = 0; i < image->numCols; i++) {
    221112            if (!mask->data.U8[j][i]) {
    222                 image->data.F32[j][i] -= background->data.F32[j][i];
     113                image->data.F32[j][i] -= backData[j][i];
    223114                if (backSub) {
    224                     backSub->data.F32[j][i] = image->data.F32[j][i];
     115                    backSub->image->data.F32[j][i] = image->data.F32[j][i];
    225116                }
    226117            }
    … …  
    229120
    230121    psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot"));
    231     psFree (rnd);
     122    psImageClippedStatsCleanup();
    232123
    233     // it is safe to free all of these: either they are saved on config->files or are local
    234     psFree (model);
    235     psFree (background);
    236     psFree (backSub);
     124    // the pmReadout selected in this function are all view on entries in config->files
    237125    return true;
    238126}
    239 
    240 double psImageClippedStats (psImage *image, psImage *mask, psU8 maskValue, double fmin, double fmax) {
    241 
    242     double value;
    243     int nx = image->numCols;
    244     int ny = image->numRows;
    245    
    246     if (nx*ny <= 0) return 0.0;
    247 
    248     int Nsubset = PS_MIN (MAX_SAMPLE_PIXELS, nx*ny);
    249     int Npixels = nx*ny;
    250 
    251     psVector *values = psVectorAlloc (Nsubset, PS_TYPE_F32);
    252 
    253     float min = values->data.F32[0];
    254     float max = values->data.F32[0];
    255 
    256     int n = 0;
    257     for (int i = 0; i < Nsubset; i++) {
    258         double frnd = psRandomUniform (rnd);
    259         int pixel = Npixels * frnd;
    260         int ix = pixel % nx;
    261         int iy = pixel / nx;
    262 
    263         if (mask->data.U8[iy][ix] & maskValue) continue;
    264 
    265         value = image->data.F32[iy][ix];
    266         min = PS_MIN (value, min);
    267         max = PS_MIN (value, max);
    268         values->data.F32[n] = value;
    269         n++;
    270     }
    271     values->n = n;
    272 
    273     int imin = fmin * n;
    274     int imax = fmax * n;
    275     int npts = imax - imin + 1;
    276 
    277     fsort (values->data.F32, n);
    278 
    279     value = 0;
    280     for (int i = imin; (i <= imax) && (i < n); i++) {
    281         value += values->data.F32[i];
    282     }
    283     value = value / npts;
    284    
    285     // XXX correct for selection bias??
    286 
    287     psFree (values);
    288     return value;
    289 }
    290 
    291 void fsort (float *value, int N) {
    292 
    293     int l,j,ir,i;
    294     float temp;
    295  
    296     if (N < 2) return;
    297     l = N >> 1;
    298     ir = N - 1;
    299     for (;;) {
    300         if (l > 0) {
    301             temp = value[--l];
    302         }
    303         else {
    304             temp = value[ir];
    305             value[ir] = value[0];
    306             if (--ir == 0) {
    307                 value[0] = temp;
    308                 return;
    309             }
    310         }
    311         i = l;
    312         j = (l << 1) + 1;
    313         while (j <= ir) {
    314             if (j < ir && value[j] < value[j+1]) ++j;
    315             if (temp < value[j]) {
    316                 value[i]=value[j];
    317                 j += (i=j) + 1;
    318             }
    319             else j = ir + 1;
    320         }
    321         value[i] = temp;
    322     }
    323 }
    324 
Note: See TracChangeset for help on using the changeset viewer.