IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27270


Ignore:
Timestamp:
Mar 12, 2010, 9:38:55 AM (16 years ago)
Author:
eugene
Message:

lots of work on source size / cr masking -- still not done!

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20100225/psphot/src/psphotSourceSize.c

    r26894 r27270  
    1313    float soft;
    1414    int grow;
     15    int xtest, ytest;
    1516} psphotSourceSizeOptions;
    1617
    … …  
    101102    assert (status);
    102103
     104    // XXX recipe name is not great
     105    options.xtest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.XTEST");
     106    options.ytest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.YTEST");
     107    assert (status);
     108
    103109    options.grow = psMetadataLookupS32(&status, recipe, "PSPHOT.CR.GROW"); // Growth size for CRs
    104110    if (!status || options.grow < 0) {
    … …  
    239245            continue;
    240246        }
     247        // psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
    241248    }
    242249
    243250    return true;
    244251}
     252
     253# define SIZE_SN_LIM 10
    245254
    246255bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
    … …  
    285294        }
    286295
    287         // we are basically classifying by moments; use the default if not found
     296        // we are basically classifying by moments
    288297        psAssert (source->moments, "why is this source missing moments?");
    289298        if (source->mode & noMoments) {
    … …  
    292301        }
    293302
     303        // convert to Mmaj, Mmin:
    294304        psF32 Mxx = source->moments->Mxx;
    295305        psF32 Myy = source->moments->Myy;
    … …  
    315325        float apMag = -2.5*log10(source->moments->Sum);
    316326        float dMag = source->psfMag - apMag;
    317         float nSigma = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);
    318 
    319         source->extNsigma = nSigma;
    320         source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     327
     328        // set nSigma to include both systematic and poisson error terms
     329        // XXX the 'poisson error' contribution for size is probably wrong...
     330        float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);
     331        float nSigmaMXX = (Mxx - psfClump->X) / hypot(psfClump->dX, psfClump->X*psfClump->X*source->errMag);
     332        float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag);
     333
     334        // partially-masked sources are more likely to be mis-measured PSFs
     335        float sizeBias = 1.0;
     336        if (source->pixWeight < 0.9) {
     337            sizeBias = 3.0;
     338        }
     339
     340        float minMxx = psfClump->X - sizeBias*options->nSigmaMoments*psfClump->dX;
     341        float minMyy = psfClump->Y - sizeBias*options->nSigmaMoments*psfClump->dY;
     342
     343        // include MAG, MXX, and MYY?
     344        source->extNsigma = nSigmaMAG;
     345
     346        // notes to clarify the source size classification rules:
     347        // * a defect should be functionally equivalent to a cosmic ray
     348        // * CR & defect should have a faintess limit (min S/N)
     349        // * SAT stars should not be faint, but defects may?
    321350
    322351        // Anything within this region is a probably PSF-like object. Saturated stars may land
    323352        // in this region, but are detected elsewhere on the basis of their peak value.
    324         bool isPSF = ((fabs(nSigma) < options->nSigmaApResid) &&
    325                       (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) &&
    326                       (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY));
     353        bool isPSF = (fabs(nSigmaMAG) < options->nSigmaApResid) && (fabs(nSigmaMXX) < sizeBias*options->nSigmaMoments) && (fabs(nSigmaMYY) < sizeBias*options->nSigmaMoments);
    327354        if (isPSF) {
    328           psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
    329                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
    330                   options->nSigmaApResid,options->nSigmaMoments);
     355          psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
     356                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
     357                  options->nSigmaApResid,sizeBias*options->nSigmaMoments);
     358          source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    331359          Npsf ++;
    332360          continue;
    … …  
    336364        // Defects may also be marked as SATSTAR -- XXX deactivate this flag?
    337365        // XXX this rule is not great
    338         if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) {
    339 
    340           psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g CR\t%g %g\n",
    341                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
    342                   options->nSigmaApResid,options->nSigmaMoments);
     366        // XXX only accept brightish detections as CRs
     367        // (nSigmaMAG < -options->nSigmaApResid) ||
     368        bool isCR = isCR = (source->errMag < 1.0 / SIZE_SN_LIM) && ((Mxx < minMxx) || (Myy < minMyy));
     369        if (isCR) {
     370            psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g %f\t%g %g  %g %g  %g %g\t%g %g\t%g CR\t%g %g\n",
     371                    source->peak->xf,source->peak->yf,source->pixWeight,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
     372                    options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    343373            source->mode |= PM_SOURCE_MODE_DEFECT;
     374            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_CR_CANDIDATE;
    344375            Ncr ++;
    345376            continue;
    … …  
    349380        // just large saturated regions.
    350381        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    351           psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g SAT\t%g %g\n",
    352                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
    353                   options->nSigmaApResid,options->nSigmaMoments);
    354          
     382            psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g SAT\t%g %g\n",
     383                    source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
     384                    options->nSigmaApResid,sizeBias*options->nSigmaMoments);
     385            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    355386            Nsat ++;
    356387            continue;
    … …  
    358389
    359390        // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)?
    360         bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y));
     391        bool isEXT = (nSigmaMAG > options->nSigmaApResid) || (Mxx > minMxx) || (Myy > minMyy);
    361392        if (isEXT) {
    362           psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Ext\t%g %g\n",
    363                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
    364                   options->nSigmaApResid,options->nSigmaMoments);
     393          psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Ext\t%g %g\n",
     394                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
     395                  options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    365396         
    366397            source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
     398            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    367399            Next ++;
    368400            continue;
    369401        }
    370         psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Unk\t%g %g\n",
    371                 source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
    372                 options->nSigmaApResid,options->nSigmaMoments);
     402        psTrace("psphotSourceClassRegion.MISS",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Unk\t%g %g\n",
     403                source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
     404                options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    373405       
    374         psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma);
     406        psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigmaMAG);
    375407        Nmiss ++;
    376408    }
    … …  
    385417bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
    386418
     419    psTimerStart ("psphot.cr");
     420
     421    int nMasked = 0;
    387422    for (int i = 0; i < sources->n; i++) {
    388423        pmSource *source = sources->data[i];
    … …  
    393428            continue;
    394429        }
     430
     431        // only check candidates marked above
     432        if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_CR_CANDIDATE)) {
     433            psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
     434            continue;
     435        }
     436
     437        // skip unless this source is thought to be a cosmic ray.  flag the detection and mask the pixels
     438        // XXX this may be degenerate with the above test
     439        if (!(source->mode & PM_SOURCE_MODE_DEFECT)) continue;
    395440
    396441        // Integer position of peak
    … …  
    402447            yPeak < 1 || yPeak > source->pixels->numRows - 2) {
    403448            psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
    404             //      psTrace("psphot.czw", 2, "Not calculating crNsigma due to edge\n");
    405449            continue;
    406450        }
    407451       
    408         // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
    409         if (source->mode & PM_SOURCE_MODE_DEFECT) {
    410             // XXX this is running slowly and is too agressive, but it more-or-less works
    411             // psphotMaskCosmicRayCZW(readout, source, options->crMask);
    412            
    413             // XXX these are the old versions which we are not using
    414             // psphotMaskCosmicRay (readout->mask, source, maskVal, crMask);
    415             // psphotMaskCosmicRay_Old (source, maskVal, crMask);
    416         }
     452        // XXX for testing, only CRMASK a single source:
     453        if (options->xtest && (fabs(source->peak->xf - options->xtest) > 5)) continue;
     454        if (options->ytest && (fabs(source->peak->yf - options->ytest) > 5)) continue;
     455
     456        // replace object in image
     457        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
     458            pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
     459        }
     460
     461        // XXX this is running slowly and is too agressive, but it more-or-less works
     462        psTrace("psphot", 6, "mask cosmic ray at %f, %f\n", source->peak->xf, source->peak->yf);
     463        psphotMaskCosmicRay(readout, source, options->crMask);
     464        nMasked ++;
     465
     466        // re-subtract the object, leave local sky
     467        pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
    417468    }
    418469   
    … …  
    430481    }
    431482
     483    psLogMsg ("psphot.cr", PS_LOG_INFO, "mask CR: %d masked in %f sec\n", nMasked, psTimerMark ("psphot.cr"));
     484
    432485    // XXX test : save the mask image
    433     if (0) {
     486    if (1) {
    434487        psphotSaveImage (NULL, readout->mask,   "mask.fits");
    435488    }
     489
    436490    return true;
    437491}
     492
     493# define LIMIT_XRANGE(X, IMAGE) { X = PS_MIN(PS_MAX(0, X), IMAGE->numCols); }
     494# define LIMIT_YRANGE(Y, IMAGE) { Y = PS_MIN(PS_MAX(0, Y), IMAGE->numRows); }
    438495
    439496// Comments by CZW 20091209 : Mechanics of how to identify CR pixels taken from "Cosmic-Ray
    … …  
    448505    pmFootprint *footprint = peak->footprint;
    449506
    450     int xm = footprint->bbox.x0;
    451     int xM = footprint->bbox.x1;
    452     int ym = footprint->bbox.y0;
    453     int yM = footprint->bbox.y1;
    454 
    455     if (xm < 0) { xm = 0; }
    456     if (ym < 0) { ym = 0; }
    457     if (xM > mask->numCols) { xM = mask->numCols; }
    458     if (yM > mask->numRows) { yM = mask->numRows; }
    459     int dx = xM - xm;
    460     int dy = yM - ym;
    461 
    462507    // Bounding boxes are inclusive of final pixel
    463     // XXX ensure dx,dy do not become too large here
    464     dx++;
    465     dy++;
     508    int xs = footprint->bbox.x0;
     509    int xe = footprint->bbox.x1 + 1;
     510    int ys = footprint->bbox.y0;
     511    int ye = footprint->bbox.y1 + 1;
     512
     513    LIMIT_XRANGE(xs, mask);
     514    LIMIT_XRANGE(xe, mask);
     515    LIMIT_YRANGE(ys, mask);
     516    LIMIT_YRANGE(ye, mask);
     517
     518    int dx = xe - xs;
     519    int dy = ye - ys;
    466520
    467521    psImage *image= readout->image;
    468522    psImage *variance = readout->variance;
    469523     
    470     int binning = 1;
    471     float sigma_thresh = 2.0;
    472     int iteration = 0;
     524    int binning = 2;
     525    float sigma_thresh = 3.0;
    473526    int max_iter = 5;
    474527    float noise_factor = 5.0 / 4.0;  // Intrinsic to the Laplacian making noise spikes spikier.
    … …  
    476529    // Temporary images.
    477530    psImage *mypix  = psImageAlloc(dx,dy,image->type.type);
     531    psImage *myfix  = psImageAlloc(dx,dy,image->type.type);
    478532    psImage *myvar  = psImageAlloc(dx,dy,image->type.type);
    479533    psImage *binned = psImageAlloc(dx * binning,dy * binning,image->type.type);
    … …  
    482536    psImage *mymask = psImageAlloc(dx,dy,PS_TYPE_IMAGE_MASK);
    483537     
    484     int x,y;
    485538    // Load my copy of things.
    486     for (x = 0; x < dx; x++) {
    487         for (y = 0; y < dy; y++) {
    488             psImageSet(mypix,x,y,psImageGet(image,x+xm,y+ym));
    489             psImageSet(myvar,x,y,psImageGet(variance,x+xm,y+ym));
     539    for (int y = 0; y < dy; y++) {
     540        for (int x = 0; x < dx; x++) {
     541            mypix->data.F32[y][x] = image->data.F32[y+ys][x+xs];
     542            myvar->data.F32[y][x] = variance->data.F32[y+ys][x+xs];
    490543            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0x00;
    491544        }
    … …  
    495548        pmSpan *sp = footprint->spans->data[i];
    496549        for (int j = sp->x0; j <= sp->x1; j++) {
    497             y = sp->y - ym;
    498             x = j - xm;
     550            int y = sp->y - ys;
     551            int x = j - xs;
    499552            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x01;
    500553        }
    501554    }
    502555
    503     int CRpix_count = 0;     
    504     do {
    505         CRpix_count = 0;
    506         // Zero out my temp images.
    507         for (x = 0; x < binning * dx; x++) {
    508             for (y = 0; y < binning * dy; y++) {
    509                 psImageSet(binned,x,y,0.0);
    510                 psImageSet(conved,x,y,0.0);
    511                 psImageSet(edges,x/binning,y/ binning,0.0);
     556    int nCRpix = 1; // force at least one pass...
     557    for (int iteration = 0; (iteration < max_iter) && (nCRpix > 0); iteration++) {
     558        nCRpix = 0;
     559        psImageInit (binned, 0.0);
     560        psImageInit (conved, 0.0);
     561        psImageInit (edges, 0.0);
     562
     563        // Make subsampled image. Maybe this should be called "unbinned" or something
     564        for (int y = 0; y < binning * dy; y++) {
     565            int yraw = y / binning;
     566            for (int x = 0; x < binning * dx; x++) {
     567                int xraw = x / binning;
     568                binned->data.F32[y][x] = mypix->data.F32[yraw][xraw];
    512569            }
    513         }     
    514         // Make subsampled image. Maybe this should be called "unbinned" or something
    515         for (x = 0; x < binning * dx; x++) {
    516             for (y = 0; y < binning * dy; y++) {
    517                 psImageSet(binned,x,y,psImageGet(mypix,x / binning,y / binning));
     570        }
     571
     572        // Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero
     573        for (int y = 1; y < binning * dy - 1; y++) {
     574            for (int x = 1; x < binning * dx - 1; x++) {
     575                float value = binned->data.F32[y][x] - 0.25 *
     576                    (binned->data.F32[y+0][x-1] + binned->data.F32[y+0][x+1] +
     577                     binned->data.F32[y-1][x+0] + binned->data.F32[y+1][x+0]);
     578                value = PS_MAX(0.0, value);
     579
     580                conved->data.F32[y][x] = value;
    518581            }
    519582        }
    520         // Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero
    521         for (x = 1; x < dx - 1; x++) {
    522             for (y = 1; y < dy - 1; y++) {
    523                 psImageSet(conved,x,y,psImageGet(binned,x,y) - 0.25 *
    524                            (psImageGet(binned,x-1,y) + psImageGet(binned,x+1,y) +
    525                             psImageGet(binned,x,y-1) + psImageGet(binned,x,y+1)));
    526                 if (psImageGet(conved,x,y) < 0.0) {
    527                     psImageSet(conved,x,y,0.0);
     583
     584        // Create an edge map by rebinning
     585        for (int y = 0; y < binning * dy; y++) {
     586            int yraw = y / binning;
     587            for (int x = 0; x < binning * dx; x++) {
     588                int xraw = x / binning;
     589                edges->data.F32[yraw][xraw] += conved->data.F32[y][x];
     590            }
     591        }
     592
     593        // Modify my mask if we're above the significance threshold
     594        for (int y = 0; y < dy; y++) {
     595            for (int x = 0; x < dx; x++) {
     596                if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x01)) continue;
     597
     598                float noise = binning * sqrt(noise_factor * myvar->data.F32[y][x]);
     599                float value = edges->data.F32[y][x] / noise;
     600
     601                if (value > sigma_thresh ) {
     602                    mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x40;
     603                    nCRpix ++;
    528604                }
    529605            }
    530606        }
    531         // Create an edge map by rebinning
    532         for (x = 0; x < binning * dx; x++) {
    533             for (y = 0; y < binning * dy; y++) {
    534                 psImageSet(edges,x / binning, y / binning,
    535                            psImageGet(edges, x / binning, y / binning) +
    536                            psImageGet(conved,x,y));
    537             }
    538         }
    539         // Modify my mask if we're above the significance threshold
    540         for (x = 0; x < dx; x++) {
    541             for (y = 0; y < dy; y++) {
    542                 if ( psImageGet(edges,x,y) / (binning * sqrt(noise_factor * psImageGet(myvar,x,y))) > sigma_thresh ) {
    543                     if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x01) {
    544                         mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x40;
    545                         CRpix_count++;
    546                     }
    547                 }
    548             }
    549         }
     607
     608# if 1
     609        psphotSaveImage (NULL, mypix,   "crmask.pix.fits");
     610# endif
    550611
    551612        // "Repair" Masked pixels for the next round.
    552         for (x = 1; x < dx - 1; x++) {
    553             for (y = 1; y < dy - 1; y++) {
    554                 if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
    555                     psImageSet(mypix,x,y,0.25 *
    556                                (psImageGet(mypix,x-1,y) + psImageGet(mypix,x+1,y) +
    557                                 psImageGet(mypix,x,y-1) + psImageGet(mypix,x,y+1)));
    558                 }
    559             }
    560         }
    561        
    562 # if 0
    563         if ((psTraceGetLevel("psphot.czw") >= 2)&&(abs(xm - 2770) < 5)&&(abs(ym - 2581) < 5)&&(iteration == 0)) {
    564           psTrace("psphot.czw",2,"TRACEMOTRON %d %d %d %d %d\n",xm,ym,dx,dy,iteration);
    565           psphotSaveImage (NULL, mypix,   "czw.pix.fits");
    566           psphotSaveImage (NULL, myvar,   "czw.var.fits");
    567           psphotSaveImage (NULL, binned,  "czw.binn.fits");
    568           psphotSaveImage (NULL, conved,  "czw.conv.fits");
    569           psphotSaveImage (NULL, edges,   "czw.edge.fits");
    570           psphotSaveImage (NULL, mymask,  "czw.mask.fits");
    571         }
    572 # endif
    573 
    574         psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration,CRpix_count);
    575         iteration++;
    576     } while ((iteration < max_iter)&&(CRpix_count > 0));
    577 
    578     // A solitary masked pixel is likely a lie. Remove those.
    579     for (x = 0; x < dx; x++) {
    580         for (y = 0; y < dy; y++) {
    581             if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
    582                 if ((x-1 >= 0)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) {
     613        for (int y = 1; y < dy - 1; y++) {
     614            for (int x = 1; x < dx - 1; x++) {
     615                if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) {
     616                    myfix->data.F32[y][x] = mypix->data.F32[y][x];
    583617                    continue;
    584618                }
    585                 if ((y-1 >= 0)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) {
    586                     continue;
    587                 }
    588                 if ((x+1 < dx)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) {
    589                     continue;
    590                 }
    591                 if ((y+1 < dy)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) {
    592                     continue;
    593                 }
    594                 mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40;
     619                myfix->data.F32[y][x] = 0.25 *
     620                    (mypix->data.F32[y+0][x-1] + mypix->data.F32[y+0][x+1] +
     621                     mypix->data.F32[y-1][x+0] + mypix->data.F32[y+1][x+0]);
    595622            }
    596623        }
    597     }
    598 
    599     // transfer temporary mask to real mask
    600     for (x = 0; x < dx; x++) {
    601         for (y = 0; y < dy; y++) {
     624       
     625        // "Repair" Masked pixels for the next round.
     626        for (int y = 1; y < dy - 1; y++) {
     627            for (int x = 1; x < dx - 1; x++) {
     628                mypix->data.F32[y][x] = myfix->data.F32[y][x];
     629            }
     630        }
     631
     632# if 1
     633        fprintf (stderr, "CRMASK %d %d %d %d %d\n", xs, ys, dx, dy, iteration);
     634        psphotSaveImage (NULL, mypix,   "crmask.fix.fits");
     635        psphotSaveImage (NULL, myvar,   "crmask.var.fits");
     636        psphotSaveImage (NULL, binned,  "crmask.binn.fits");
     637        psphotSaveImage (NULL, conved,  "crmask.conv.fits");
     638        psphotSaveImage (NULL, edges,   "crmask.edge.fits");
     639        psphotSaveImage (NULL, mymask,  "crmask.mask.fits");
     640# endif
     641        psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration, nCRpix);
     642    }
     643
     644    // A solitary masked pixel is likely a lie. Remove those
     645    // XXX can't we use nCRpix == 1 to test for these?
     646    for (int x = 0; x < dx; x++) {
     647        for (int y = 0; y < dy; y++) {
     648            if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) continue;
     649            if ((x-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) {
     650                continue;
     651            }
     652            if ((y-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) {
     653                continue;
     654            }
     655            if ((x+1 < dx) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) {
     656                continue;
     657            }
     658            if ((y+1 < dy) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) {
     659                continue;
     660            }
     661            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40;
     662        }
     663    }
     664
     665    // transfer temporary mask to real mask & count masked pixels
     666    nCRpix = 0;
     667    for (int x = 0; x < dx; x++) {
     668        for (int y = 0; y < dy; y++) {
    602669            if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
    603                 mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ym+mask->row0][x+xm+mask->col0] |= maskVal;
     670                mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ys+mask->row0][x+xs+mask->col0] |= maskVal;
     671                nCRpix ++;
    604672            }
    605673        }
    … …  
    607675     
    608676    // XXX if we decide this REALLY is a cosmic ray, set the CR_LIMIT bit
    609     source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     677    if (nCRpix > 1) {
     678        source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     679    }
    610680
    611681    psFree(mypix);
     682    psFree(myfix);
    612683    psFree(myvar);
    613684    psFree(binned);
Note: See TracChangeset for help on using the changeset viewer.