IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 18868


Ignore:
Timestamp:
Aug 1, 2008, 3:29:28 PM (18 years ago)
Author:
Paul Price
Message:

Grow the CR mask.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotSourceSize.c

    r18835 r18868  
    1010// deviation from the psf model at the r = FWHM/2 position
    1111
    12 bool psphotSourceSize (pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, long first) {
     12bool psphotSourceSize (pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, long first)
     13{
    1314
    1415    bool status;
    … …  
    2930    assert (status);
    3031
     32    int grow = psMetadataLookupS32(&status, recipe, "PSPHOT.CR.GROW"); // Growth size for CRs
     33    if (!status || grow < 0) {
     34        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSPHOT.CR.GROW is not positive.");
     35        return false;
     36    }
     37
    3138    // loop over all source
    3239    for (int i = first; i < sources->n; i++) {
    33         pmSource *source = sources->data[i];
    34 
    35         // skip source if it was already measured
    36         if (isfinite(source->crNsigma)) continue;
    37 
    38         // source must have been subtracted
    39         if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue;
    40 
    41         psF32 **resid  = source->pixels->data.F32;
    42         psF32 **weight = source->weight->data.F32;
    43         psU8 **mask    = source->maskObj->data.U8;
    44 
    45         // check for extendedness: measure the delta flux significance at the 1 sigma contour
    46         source->extNsigma = psphotModelContour (source->pixels, source->weight, source->maskObj, source->modelPSF, 1.0);
    47 
    48         // XXX prevent a source from being both CR and EXT?
    49         if (source->extNsigma > EXT_NSIGMA_LIMIT) {
    50           source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
    51         }
    52 
    53         int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
    54         int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
    55 
    56         // XXX for now, skip sources which are too close to a boundary
    57         // XXX raise a flag?
    58         if (xPeak < 1) continue;
    59         if (xPeak > source->pixels->numCols - 2) continue;
    60         if (yPeak < 1) continue;
    61         if (yPeak > source->pixels->numRows - 2) continue;
    62 
    63         // XXX for now, just skip any sources with masked pixels
    64         // XXX raise a flag?
    65         bool keep = true;
    66         for (int iy = -1; (iy <= +1) && keep; iy++) {
    67             for (int ix = -1; (ix <= +1) && keep; ix++) {
    68                 if (mask[yPeak+iy][xPeak+ix]) { keep &= false; }
    69             }
    70         }
    71         if (!keep) continue;
    72 
    73         // XXX need to deal with edge peaks... and mask
    74         float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1];
    75         float cY = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak+0] - resid[yPeak+1][xPeak+0];
    76         float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1];
    77         float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1];
    78 
    79         float dcX = 4*weight[yPeak][xPeak] + weight[yPeak+0][xPeak-1] + weight[yPeak+0][xPeak+1];
    80         float dcY = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak+0] + weight[yPeak+1][xPeak+0];
    81         float dcL = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak-1] + weight[yPeak+1][xPeak+1];
    82         float dcR = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak-1] + weight[yPeak-1][xPeak+1];
    83 
    84         float nX = cX / sqrt(dcX);
    85         float nY = cY / sqrt(dcY);
    86         float nL = cL / sqrt(dcL);
    87         float nR = cR / sqrt(dcR);
    88 
    89         // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
    90         // Ndof = 4 ? (four measurements, no free parameters)
    91         // XXX this value is going to be biased low because of systematic errors.
    92         // we need to calibrate it somehow
    93         // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq);
    94 
    95         // not strictly accurate: overcounts the chisq contribution from the center pixel (by factor of 4)
    96         source->psfChisq = PS_SQR (nX) + PS_SQR (nY) + PS_SQR (nX) + PS_SQR (nR);
    97 
    98         float fCR = 0.0;
    99         float fEXT = 0.0;
    100         int nCR = 0;
    101         int nEXT = 0;
    102         if (nX > 0.0) {
    103             fCR += nX;
    104             nCR ++;
    105         } else {
    106             fEXT += nX;
    107             nEXT ++;
    108         }
    109         if (nY > 0.0) {
    110             fCR += nY;
    111             nCR ++;
    112         } else {
    113             fEXT += nY;
    114             nEXT ++;
    115         }
    116         if (nL > 0.0) {
    117             fCR += nL;
    118             nCR ++;
    119         } else {
    120             fEXT += nL;
    121             nEXT ++;
    122         }
    123         if (nR > 0.0) {
    124             fCR += nR;
    125             nCR ++;
    126         } else {
    127             fEXT += nR;
    128             nEXT ++;
    129         }
    130         source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
    131         // NOTE: abs needed to make the Nsigma value positive
    132 
    133         if (!isfinite(source->crNsigma)) {
    134           fprintf (stderr, ".");
    135           source->crNsigma = -6.0;
    136         }
    137 
    138         // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
    139         if (source->crNsigma > CR_NSIGMA_LIMIT) {
    140             source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    141             pmPeak *peak = source->peak;
    142             psAssert (peak, "NULL peak");
    143 
    144             // replace the source flux
    145             pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    146             source->mode &= ~PM_SOURCE_MODE_SUBTRACTED;
    147 
    148             psImage *mask   = source->maskView;
    149             psImage *pixels = source->pixels;
    150             psImage *weight = source->weight;
    151 
    152             # define SN_LIMIT 5.0
    153 
    154             int xo = peak->x - pixels->col0;
    155             int yo = peak->y - pixels->row0;
    156 
    157             // mark the pixels in this row to the left, then the right
    158             for (int ix = xo; ix >= 0; ix--) {
    159                 float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]);
    160                 if (SN > SN_LIMIT) {
    161                     mask->data.U8[yo][ix] |= crMask;
    162                 }
    163             }
    164             for (int ix = xo + 1; ix < pixels->numCols; ix++) {
    165                 float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]);
    166                 if (SN > SN_LIMIT) {
    167                     mask->data.U8[yo][ix] |= crMask;
    168                 }
    169             }
    170 
    171             // for each of the neighboring rows, mark the high pixels if they have a marked neighbor
    172             // first go up:
    173             for (int iy = yo; iy >= 0; iy--) {
    174                 // mark the pixels in this row to the left, then the right
    175                 for (int ix = 0; ix < pixels->numCols; ix++) {
    176                     float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
    177                     if (SN < SN_LIMIT) continue;
    178                    
    179                     bool valid = false;
    180                     valid |= (mask->data.U8[iy+1][ix] & crMask);
    181                     valid |= (ix > 0) ? (mask->data.U8[iy+1][ix-1] & crMask) : 0;
    182                     valid |= (ix <= mask->numCols) ? (mask->data.U8[iy+1][ix+1] & crMask) : 0;
    183 
    184                     if (!valid) continue;
    185                     mask->data.U8[iy][ix] |= crMask;
    186                 }
    187             }
    188             // next go down:
    189             for (int iy = yo+1; iy < pixels->numRows; iy++) {
    190                 // mark the pixels in this row to the left, then the right
    191                 for (int ix = 0; ix < pixels->numCols; ix++) {
    192                     float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
    193                     if (SN < SN_LIMIT) continue;
    194                    
    195                     bool valid = false;
    196                     valid |= (mask->data.U8[iy-1][ix] & crMask);
    197                     valid |= (ix > 0) ? (mask->data.U8[iy-1][ix-1] & crMask) : 0;
    198                     valid |= (ix <= mask->numCols) ? (mask->data.U8[iy-1][ix+1] & crMask) : 0;
    199 
    200                     if (!valid) continue;
    201                     mask->data.U8[iy][ix] |= crMask;
    202                 }
    203             }
    204         }
     40        pmSource *source = sources->data[i];
     41
     42        // skip source if it was already measured
     43        if (isfinite(source->crNsigma)) continue;
     44
     45        // source must have been subtracted
     46        if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue;
     47
     48        psF32 **resid  = source->pixels->data.F32;
     49        psF32 **weight = source->weight->data.F32;
     50        psU8 **mask    = source->maskObj->data.U8;
     51
     52        // check for extendedness: measure the delta flux significance at the 1 sigma contour
     53        source->extNsigma = psphotModelContour (source->pixels, source->weight, source->maskObj, source->modelPSF, 1.0);
     54
     55        // XXX prevent a source from being both CR and EXT?
     56        if (source->extNsigma > EXT_NSIGMA_LIMIT) {
     57          source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
     58        }
     59
     60        int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
     61        int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
     62
     63        // XXX for now, skip sources which are too close to a boundary
     64        // XXX raise a flag?
     65        if (xPeak < 1) continue;
     66        if (xPeak > source->pixels->numCols - 2) continue;
     67        if (yPeak < 1) continue;
     68        if (yPeak > source->pixels->numRows - 2) continue;
     69
     70        // XXX for now, just skip any sources with masked pixels
     71        // XXX raise a flag?
     72        bool keep = true;
     73        for (int iy = -1; (iy <= +1) && keep; iy++) {
     74            for (int ix = -1; (ix <= +1) && keep; ix++) {
     75                if (mask[yPeak+iy][xPeak+ix]) { keep &= false; }
     76            }
     77        }
     78        if (!keep) continue;
     79
     80        // XXX need to deal with edge peaks... and mask
     81        float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1];
     82        float cY = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak+0] - resid[yPeak+1][xPeak+0];
     83        float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1];
     84        float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1];
     85
     86        float dcX = 4*weight[yPeak][xPeak] + weight[yPeak+0][xPeak-1] + weight[yPeak+0][xPeak+1];
     87        float dcY = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak+0] + weight[yPeak+1][xPeak+0];
     88        float dcL = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak-1] + weight[yPeak+1][xPeak+1];
     89        float dcR = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak-1] + weight[yPeak-1][xPeak+1];
     90
     91        float nX = cX / sqrt(dcX);
     92        float nY = cY / sqrt(dcY);
     93        float nL = cL / sqrt(dcL);
     94        float nR = cR / sqrt(dcR);
     95
     96        // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
     97        // Ndof = 4 ? (four measurements, no free parameters)
     98        // XXX this value is going to be biased low because of systematic errors.
     99        // we need to calibrate it somehow
     100        // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq);
     101
     102        // not strictly accurate: overcounts the chisq contribution from the center pixel (by factor of 4)
     103        source->psfChisq = PS_SQR (nX) + PS_SQR (nY) + PS_SQR (nX) + PS_SQR (nR);
     104
     105        float fCR = 0.0;
     106        float fEXT = 0.0;
     107        int nCR = 0;
     108        int nEXT = 0;
     109        if (nX > 0.0) {
     110            fCR += nX;
     111            nCR ++;
     112        } else {
     113            fEXT += nX;
     114            nEXT ++;
     115        }
     116        if (nY > 0.0) {
     117            fCR += nY;
     118            nCR ++;
     119        } else {
     120            fEXT += nY;
     121            nEXT ++;
     122        }
     123        if (nL > 0.0) {
     124            fCR += nL;
     125            nCR ++;
     126        } else {
     127            fEXT += nL;
     128            nEXT ++;
     129        }
     130        if (nR > 0.0) {
     131            fCR += nR;
     132            nCR ++;
     133        } else {
     134            fEXT += nR;
     135            nEXT ++;
     136        }
     137        source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
     138        // NOTE: abs needed to make the Nsigma value positive
     139
     140        if (!isfinite(source->crNsigma)) {
     141          fprintf (stderr, ".");
     142          source->crNsigma = -6.0;
     143        }
     144
     145        // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
     146        if (source->crNsigma > CR_NSIGMA_LIMIT) {
     147            source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     148            pmPeak *peak = source->peak;
     149            psAssert (peak, "NULL peak");
     150
     151            // replace the source flux
     152            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     153            source->mode &= ~PM_SOURCE_MODE_SUBTRACTED;
     154
     155            psImage *mask   = source->maskView;
     156            psImage *pixels = source->pixels;
     157            psImage *weight = source->weight;
     158
     159            # define SN_LIMIT 5.0
     160
     161            int xo = peak->x - pixels->col0;
     162            int yo = peak->y - pixels->row0;
     163
     164            // mark the pixels in this row to the left, then the right
     165            for (int ix = xo; ix >= 0; ix--) {
     166                float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]);
     167                if (SN > SN_LIMIT) {
     168                    mask->data.U8[yo][ix] |= crMask;
     169                }
     170            }
     171            for (int ix = xo + 1; ix < pixels->numCols; ix++) {
     172                float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]);
     173                if (SN > SN_LIMIT) {
     174                    mask->data.U8[yo][ix] |= crMask;
     175                }
     176            }
     177
     178            // for each of the neighboring rows, mark the high pixels if they have a marked neighbor
     179            // first go up:
     180            for (int iy = yo; iy >= 0; iy--) {
     181                // mark the pixels in this row to the left, then the right
     182                for (int ix = 0; ix < pixels->numCols; ix++) {
     183                    float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
     184                    if (SN < SN_LIMIT) continue;
     185
     186                    bool valid = false;
     187                    valid |= (mask->data.U8[iy+1][ix] & crMask);
     188                    valid |= (ix > 0) ? (mask->data.U8[iy+1][ix-1] & crMask) : 0;
     189                    valid |= (ix <= mask->numCols) ? (mask->data.U8[iy+1][ix+1] & crMask) : 0;
     190
     191                    if (!valid) continue;
     192                    mask->data.U8[iy][ix] |= crMask;
     193                }
     194            }
     195            // next go down:
     196            for (int iy = yo+1; iy < pixels->numRows; iy++) {
     197                // mark the pixels in this row to the left, then the right
     198                for (int ix = 0; ix < pixels->numCols; ix++) {
     199                    float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
     200                    if (SN < SN_LIMIT) continue;
     201
     202                    bool valid = false;
     203                    valid |= (mask->data.U8[iy-1][ix] & crMask);
     204                    valid |= (ix > 0) ? (mask->data.U8[iy-1][ix-1] & crMask) : 0;
     205                    valid |= (ix <= mask->numCols) ? (mask->data.U8[iy-1][ix+1] & crMask) : 0;
     206
     207                    if (!valid) continue;
     208                    mask->data.U8[iy][ix] |= crMask;
     209                }
     210            }
     211        }
    205212    }
    206213
    207     psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n - first, psTimerMark ("psphot"));
     214    if (grow > 0 && !psImageConvolveMask(readout->mask, readout->mask, crMask, crMask,
     215                                         -grow, grow, -grow, grow)) {
     216        psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask");
     217        return false;
     218    }
     219
     220    psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n",
     221              sources->n - first, psTimerMark ("psphot"));
    208222
    209223    return true;
    … …  
    220234// any are significantly negative, some may be significantly positive : CR
    221235// any are significantly positive, none may be significantly positive : EXT
    222        
     236
    223237/* Nn  Np  No
    224238   4   0   0   CR_1
    … …  
    240254
    241255/* Alternatively, write a f(CR) = Sum(nX,etc if >0) */
    242          
     256
    243257
    244258/* I can write the formal probability that the 4 measurements are consistent with a PSF
    … …  
    273287
    274288    for (int x = xMin; x <= xMax; x++) {
    275         float A = PS_SQR (1.0 / PAR[PM_PAR_SYY]);
    276         float B = x * PAR[PM_PAR_SXY];
    277         float C = PS_SQR (x / PAR[PM_PAR_SXX]) - Ro;
    278 
    279         float T = PS_SQR(B) - 4*A*C;
    280         if (T < 0.0) continue;
    281    
    282         float yP = (-B + sqrt (T)) / (2.0 * A);
    283         float yM = (-B - sqrt (T)) / (2.0 * A);
    284 
    285         int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
    286         int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
    287         int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
    288 
    289         if (xPix < 0) continue;
    290         if (xPix >= image->numCols) continue;
    291 
    292         if ((yPixM >= 0) && (yPixM < image->numRows)) {
    293             if (!mask || !mask->data.U8[yPixM][xPix]) {
    294                 float dSigma = image->data.F32[yPixM][xPix] / sqrt (weight->data.F32[yPixM][xPix]);
    295                 nSigma += dSigma;
    296                 nPts ++;
    297             }
    298         }
    299        
    300         if (yPixM == yPixP) continue;
    301 
    302         if ((yPixP >= 0) && (yPixP < image->numRows)) {
    303             if (!mask || !mask->data.U8[yPixP][xPix]) {
    304                 float dSigma = image->data.F32[yPixP][xPix] / sqrt (weight->data.F32[yPixP][xPix]);
    305                 nSigma += dSigma;
    306                 nPts ++;
    307             }
    308         }
    309     }   
     289        float A = PS_SQR (1.0 / PAR[PM_PAR_SYY]);
     290        float B = x * PAR[PM_PAR_SXY];
     291        float C = PS_SQR (x / PAR[PM_PAR_SXX]) - Ro;
     292
     293        float T = PS_SQR(B) - 4*A*C;
     294        if (T < 0.0) continue;
     295
     296        float yP = (-B + sqrt (T)) / (2.0 * A);
     297        float yM = (-B - sqrt (T)) / (2.0 * A);
     298
     299        int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
     300        int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
     301        int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
     302
     303        if (xPix < 0) continue;
     304        if (xPix >= image->numCols) continue;
     305
     306        if ((yPixM >= 0) && (yPixM < image->numRows)) {
     307            if (!mask || !mask->data.U8[yPixM][xPix]) {
     308                float dSigma = image->data.F32[yPixM][xPix] / sqrt (weight->data.F32[yPixM][xPix]);
     309                nSigma += dSigma;
     310                nPts ++;
     311            }
     312        }
     313
     314        if (yPixM == yPixP) continue;
     315
     316        if ((yPixP >= 0) && (yPixP < image->numRows)) {
     317            if (!mask || !mask->data.U8[yPixP][xPix]) {
     318                float dSigma = image->data.F32[yPixP][xPix] / sqrt (weight->data.F32[yPixP][xPix]);
     319                nSigma += dSigma;
     320                nPts ++;
     321            }
     322        }
     323    }
    310324    nSigma /= nPts;
    311325    return nSigma;
Note: See TracChangeset for help on using the changeset viewer.