IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26618


Ignore:
Timestamp:
Jan 15, 2010, 3:55:36 PM (16 years ago)
Author:
eugene
Message:

revert to older (non czw cr-masking) version

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotSourceSize.c

    r26617 r26618  
    2222bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
    2323
    24 bool psphotMaskCosmicRayCZW (pmReadout *readout, pmSource *source, psImageMaskType maskVal,psphotSourceSizeOptions *options);
    2524// we need to call this function after sources have been fitted to the PSF model and
    2625// subtracted.  To determine the CR-nature, this function examines the 9 pixels in the 3x3
    … …  
    7170
    7271    // We are using the value psfMag - 2.5*log10(moment.Sum) as a measure of the extendedness
    73     // of an object.  We need to model this distribution for the PSF stars before we can test
     72    // of and object.  We need to model this distribution for the PSF stars before we can test
    7473    // the significance for a specific object
    7574    // XXX move this to the code that generates the PSF?
    … …  
    104103    pmFootprint *footprint = peak->footprint;
    105104    if (!footprint) {
    106       psTrace("psphot.czw",2,"Using isophot CR mask code.");
    107      
    108105        // if we have not footprint, use the old code to mask by isophot
    109106        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    … …  
    112109
    113110    if (!footprint->spans) {
    114       psTrace("psphot.czw",2,"Using isophot CR mask code.");
    115      
    116111        // if we have no footprint, use the old code to mask by isophot
    117112        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    118113        return true;
    119114    }
    120     psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    121 /*     // mask all of the pixels covered by the spans of the footprint */
    122 /*     for (int j = 1; j < footprint->spans->n; j++) { */
    123 /*         pmSpan *span1 = footprint->spans->data[j]; */
    124 
    125 /*         int iy = span1->y; */
    126 /*         int xs = span1->x0; */
    127 /*         int xe = span1->x1; */
    128 
    129 /*         for (int ix = xs; ix < xe; ix++) { */
    130 /*             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; */
    131 /*         } */
    132 /*     } */
     115
     116    // mask all of the pixels covered by the spans of the footprint
     117    for (int j = 1; j < footprint->spans->n; j++) {
     118        pmSpan *span1 = footprint->spans->data[j];
     119
     120        int iy = span1->y;
     121        int xs = span1->x0;
     122        int xe = span1->x1;
     123
     124        for (int ix = xs; ix < xe; ix++) {
     125            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     126        }
     127    }
    133128    return true;
    134129}
    135 
    136  
    137 
    138130
    139131bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
    … …  
    240232        float dMag = source->psfMag - apMag;
    241233
    242         psVectorAppend (Ap, 100, dMag);
    243         psVectorAppend (ApErr, 100, source->errMag);
     234        psVectorAppend (Ap, dMag);
     235        psVectorAppend (ApErr, source->errMag);
    244236    }
    245237
    … …  
    255247    options->ApResid = stats->robustMedian;
    256248    options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05);
     249    // XXX this is quite arbitrary...
     250    if (!isfinite(options->ApSysErr)) options->ApSysErr = 0.01;
    257251    psLogMsg ("psphot", PS_LOG_DETAIL, "psf - Sum: %f +/- %f\n", options->ApResid, options->ApSysErr);
    258252
    … …  
    381375        // Anything within this region is a probably PSF-like object. Saturated stars may land
    382376        // in this region, but are detected elsewhere on the basis of their peak value.
    383         bool isPSF = ((fabs(nSigma) < options->nSigmaApResid) &&
    384                       (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) &&
    385                       (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY));
     377        bool isPSF = (fabs(nSigma) < options->nSigmaApResid) && (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) && (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY);
    386378        if (isPSF) {
    387           psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
    388                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
    389                   options->nSigmaApResid,options->nSigmaMoments);
    390           Npsf ++;
    391           continue;
     379            Npsf ++;
     380            continue;
    392381        }
    393382
    … …  
    396385        // XXX this rule is not great
    397386        if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) {
    398 
    399           psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g CR\t%g %g\n",
    400                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
    401                   options->nSigmaApResid,options->nSigmaMoments);
    402387            source->mode |= PM_SOURCE_MODE_DEFECT;
    403388            Ncr ++;
    … …  
    408393        // just large saturated regions.
    409394        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    410           psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g SAT\t%g %g\n",
    411                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
    412                   options->nSigmaApResid,options->nSigmaMoments);
    413          
    414395            Nsat ++;
    415396            continue;
    … …  
    419400        bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y));
    420401        if (isEXT) {
    421           psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Ext\t%g %g\n",
    422                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
    423                   options->nSigmaApResid,options->nSigmaMoments);
    424          
    425402            source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
    426403            Next ++;
    427404            continue;
    428405        }
    429         psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Unk\t%g %g\n",
    430                 source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
    431                 options->nSigmaApResid,options->nSigmaMoments);
    432        
     406
    433407        psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma);
    434408        Nmiss ++;
    … …  
    521495        pmSource *source = sources->data[i];
    522496
    523 /*         // skip source if it was already measured */
    524 /*         if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) { */
    525 /*             psTrace("psphot", 7, "Not calculating source size since it has already been measured\n"); */
    526 /*          psTrace("psphot.czw", 2, "Not calculating source size since it has already been measured\n"); */
    527 /*             continue; */
    528 /*         } */
    529 
    530 /*         // source must have been subtracted */
    531 /*         if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) { */
    532 /*             source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; */
    533 /*             psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n"); */
    534 /*          psTrace("psphot.czw", 2, "Not calculating source size since source is not subtracted\n"); */
    535 /*             continue; */
    536 /*         } */
    537 
    538 //        psF32 **resid  = source->pixels->data.F32;
    539  //       psF32 **variance = source->variance->data.F32;
    540  //       psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     497        // skip source if it was already measured
     498        if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
     499            psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
     500            continue;
     501        }
     502
     503        // source must have been subtracted
     504        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
     505            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
     506            psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
     507            continue;
     508        }
     509
     510        psF32 **resid  = source->pixels->data.F32;
     511        psF32 **variance = source->variance->data.F32;
     512        psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    541513
    542514        // Integer position of peak
    … …  
    548520            yPeak < 1 || yPeak > source->pixels->numRows - 2) {
    549521            psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
    550             //      psTrace("psphot.czw", 2, "Not calculating crNsigma due to edge\n");
    551             continue;
    552         }
    553        
    554 /*         // Skip sources with masked pixels.  These are mostly caught as DEFECT */
    555 /*         bool keep = true; */
    556 /*         for (int iy = -1; (iy <= +1) && keep; iy++) { */
    557 /*             for (int ix = -1; (ix <= +1) && keep; ix++) { */
    558 /*                 if (mask[yPeak+iy][xPeak+ix] & options->maskVal) { */
    559 /*                     keep = false; */
    560 /*                 } */
    561 /*             } */
    562 /*         } */
    563 /*         if (!keep) { */
    564 /*             psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n"); */
    565 /*          //      psTrace("psphot.czw", 2, "Not calculating crNsigma due to masked pixels\n"); */
    566 /*             continue; */
    567 /*         } */
    568 
    569 /*      //      psTrace("psphot.czw", 2, "Actually doing something to mask a CR. \n"); */
    570 /*         // Compare the central pixel with those on either side, for the four possible lines through it. */
    571 
    572 /*         // Soften variances (add systematic error) */
    573 /*         float softening = options->soft * PS_SQR(source->peak->flux); // Softening for variances */
    574 
    575 /*         // Across the middle: y = 0 */
    576 /*         float cX = 2*resid[yPeak][xPeak]   - resid[yPeak+0][xPeak-1]  - resid[yPeak+0][xPeak+1]; */
    577 /*         float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1]; */
    578 /*         float nX = cX / sqrtf(dcX + softening); */
    579 
    580 /*         // Up the centre: x = 0 */
    581 /*         float cY = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak+0]  - resid[yPeak+1][xPeak+0]; */
    582 /*         float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0]; */
    583 /*         float nY = cY / sqrtf(dcY + softening); */
    584 
    585 /*         // Diagonal: x = y */
    586 /*         float cL = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak-1]  - resid[yPeak+1][xPeak+1]; */
    587 /*         float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1]; */
    588 /*         float nL = cL / sqrtf(dcL + softening); */
    589 
    590 /*         // Diagonal: x = - y */
    591 /*         float cR = 2*resid[yPeak][xPeak]   - resid[yPeak+1][xPeak-1]  - resid[yPeak-1][xPeak+1]; */
    592 /*         float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1]; */
    593 /*         float nR = cR / sqrtf(dcR + softening); */
    594 
    595 /*         // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2) */
    596 /*         // Ndof = 4 ? (four measurements, no free parameters) */
    597 /*         // XXX this value is going to be biased low because of systematic errors. */
    598 /*         // we need to calibrate it somehow */
    599 /*         // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq); */
    600 
    601 /*         // not strictly accurate: overcounts the chisq contribution from the center pixel (by */
    602 /*         // factor of 4); also biases a bit low if any pixels are masked */
    603 /*         // XXX I am not sure I want to keep this value... */
    604 /*         source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR); */
    605 
    606 /*         float fCR = 0.0; */
    607 /*         int nCR = 0; */
    608 /*         if (nX > 0.0) { */
    609 /*             fCR += nX; */
    610 /*             nCR ++; */
    611 /*         } */
    612 /*         if (nY > 0.0) { */
    613 /*             fCR += nY; */
    614 /*             nCR ++; */
    615 /*         } */
    616 /*         if (nL > 0.0) { */
    617 /*             fCR += nL; */
    618 /*             nCR ++; */
    619 /*         } */
    620 /*         if (nR > 0.0) { */
    621 /*             fCR += nR; */
    622 /*             nCR ++; */
    623 /*         } */
    624 /*         source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0; */
    625 /*         source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; */
    626 
    627 /*         if (!isfinite(source->crNsigma)) { */
    628 /*             continue; */
    629 /*         } */
    630 
    631         //psImageMaskType crMask  = 0x08;
    632        
     522            continue;
     523        }
     524
     525        // Skip sources with masked pixels.  These are mostly caught as DEFECT
     526        bool keep = true;
     527        for (int iy = -1; (iy <= +1) && keep; iy++) {
     528            for (int ix = -1; (ix <= +1) && keep; ix++) {
     529                if (mask[yPeak+iy][xPeak+ix] & options->maskVal) {
     530                    keep = false;
     531                }
     532            }
     533        }
     534        if (!keep) {
     535            psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n");
     536            continue;
     537        }
     538
     539        // Compare the central pixel with those on either side, for the four possible lines through it.
     540
     541        // Soften variances (add systematic error)
     542        float softening = options->soft * PS_SQR(source->peak->flux); // Softening for variances
     543
     544        // Across the middle: y = 0
     545        float cX = 2*resid[yPeak][xPeak]   - resid[yPeak+0][xPeak-1]  - resid[yPeak+0][xPeak+1];
     546        float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1];
     547        float nX = cX / sqrtf(dcX + softening);
     548
     549        // Up the centre: x = 0
     550        float cY = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak+0]  - resid[yPeak+1][xPeak+0];
     551        float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0];
     552        float nY = cY / sqrtf(dcY + softening);
     553
     554        // Diagonal: x = y
     555        float cL = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak-1]  - resid[yPeak+1][xPeak+1];
     556        float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1];
     557        float nL = cL / sqrtf(dcL + softening);
     558
     559        // Diagonal: x = - y
     560        float cR = 2*resid[yPeak][xPeak]   - resid[yPeak+1][xPeak-1]  - resid[yPeak-1][xPeak+1];
     561        float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1];
     562        float nR = cR / sqrtf(dcR + softening);
     563
     564        // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
     565        // Ndof = 4 ? (four measurements, no free parameters)
     566        // XXX this value is going to be biased low because of systematic errors.
     567        // we need to calibrate it somehow
     568        // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq);
     569
     570        // not strictly accurate: overcounts the chisq contribution from the center pixel (by
     571        // factor of 4); also biases a bit low if any pixels are masked
     572        // XXX I am not sure I want to keep this value...
     573        source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR);
     574
     575        float fCR = 0.0;
     576        int nCR = 0;
     577        if (nX > 0.0) {
     578            fCR += nX;
     579            nCR ++;
     580        }
     581        if (nY > 0.0) {
     582            fCR += nY;
     583            nCR ++;
     584        }
     585        if (nL > 0.0) {
     586            fCR += nL;
     587            nCR ++;
     588        }
     589        if (nR > 0.0) {
     590            fCR += nR;
     591            nCR ++;
     592        }
     593        source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
     594        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     595
     596        if (!isfinite(source->crNsigma)) {
     597            continue;
     598        }
     599
    633600        // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
    634         if (source->mode & PM_SOURCE_MODE_DEFECT) {
    635           //        if (source->crNsigma > options->nSigmaCR) {
    636           source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    637           pmPeak *peak = source->peak;
    638           pmFootprint *footprint = peak->footprint;
    639           //      if ((footprint) && (footprint->spans) && !(footprint->npix < 0)) {
    640           psTrace("psphot.czw",2,"Footprint sync: %d %d\n",peak->x,peak->y);
    641           psTrace("psphot.czw",2,"Footprint info: %p %p %d %d %d\n",(void *)footprint, (void *) footprint->spans,footprint->id,peak->x,peak->y);
    642           if ((footprint) && (footprint->spans)) {
    643             psTrace("psphot.czw",2,"  Footprint good?: %p %p %d %g %g\n",(void *)footprint, (void *) footprint->spans,
    644                     footprint->id,footprint->bbox.x0,footprint->bbox.y0);
    645           }
    646           psImageMaskType maskVal = 0xbf;  // HACK
    647           maskVal = 0x80;
    648           psphotMaskCosmicRayCZW(readout, source, maskVal,options);
    649          
    650           //                psphotMaskCosmicRay(readout->mask, source, maskVal, maskVal);
    651          
    652           // XXX still testing... : psphotMaskCosmicRay (readout->mask, source, maskVal, crMask);
    653           // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask);
    654         }
    655     }
    656    
     601        if (source->crNsigma > options->nSigmaCR) {
     602            source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     603            // XXX still testing... : psphotMaskCosmicRay (readout->mask, source, maskVal, crMask);
     604            // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask);
     605        }
     606    }
     607
    657608    // now that we have masked pixels associated with CRs, we can grow the mask
    658609    if (options->grow > 0) {
    … …  
    667618        readout->mask = newMask;
    668619    }
    669 
    670     if (psTraceGetLevel("psphot.czw") >= 2) {
    671       psphotSaveImage (NULL, readout->mask,   "mask.fits");
    672     }
    673620    return true;
    674621}
    675 
    676 
    677 // Mechanics of how to identify CR pixels taken from
    678 // "Cosmic-Ray Rejection by Laplacian Edge Detection" by Pieter van Dokkum, arXiv:astro-ph/0108003 .
    679 // This does no repair or recovery of the CR pixels, it only masks them out.
    680 // My test code can be found at /data/ipp031.0/watersc1/psphot.20091209/algo_check.c
    681 bool psphotMaskCosmicRayCZW (pmReadout *readout, pmSource *source, psImageMaskType maskVal,psphotSourceSizeOptions *options) {
    682   // Get the actual images and information about the peak.
    683   psImage *mask = readout->mask;
    684   pmPeak *peak = source->peak;
    685   pmFootprint *footprint = peak->footprint;
    686 
    687   // Try to skip sources with bad footprints.
    688   if ((footprint) && (footprint->spans) && footprint->id > 0 && footprint->id < 100000) {
    689     int xm = footprint->bbox.x0;
    690     int xM = footprint->bbox.x1;
    691     int ym = footprint->bbox.y0;
    692     int yM = footprint->bbox.y1;
    693 
    694     if (xm < 0) { xm = 0; }
    695     if (ym < 0) { ym = 0; }
    696     if (xM > mask->numCols) { xM = mask->numCols; }
    697     if (yM > mask->numRows) { yM = mask->numRows; }
    698     int dx = xM - xm;
    699     int dy = yM - ym;
    700 
    701     // Further try to skip bad footprints.
    702     if ((footprint->npix > 0) && ((dx)*(dy) > 0) && (dx < 4000) && (dy < 4000) && (dx > 0) && (dy > 0)) {
    703       // Bounding boxes are inclusive of final pixel.
    704       dx++;
    705       dy++;
    706 
    707       psImage *image= readout->image;
    708       psImage *variance = readout->variance;
    709      
    710       int binning = 1;
    711       float sigma_thresh = 2.0;
    712       int iteration = 0;
    713       int max_iter = 5;
    714       float noise_factor = 5.0 / 4.0;  // Intrinsic to the Laplacian making noise spikes spikier.
    715 
    716       // Temporary images.
    717       psImage *mypix  = psImageAlloc(dx,dy,image->type.type);
    718       psImage *myvar  = psImageAlloc(dx,dy,image->type.type);
    719       psImage *binned = psImageAlloc(dx * binning,dy * binning,image->type.type);
    720       psImage *conved = psImageAlloc(dx * binning,dy * binning,image->type.type);
    721       psImage *edges  = psImageAlloc(dx,dy,image->type.type);
    722       psImage *mymask = psImageAlloc(dx,dy,PS_TYPE_IMAGE_MASK);
    723      
    724       int x,y;
    725       // Load my copy of things.
    726       for (x = 0; x < dx; x++) {
    727         for (y = 0; y < dy; y++) {
    728           psImageSet(mypix,x,y,psImageGet(image,x+xm,y+ym));
    729           psImageSet(myvar,x,y,psImageGet(variance,x+xm,y+ym));
    730           mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0x00;
    731         }
    732       }
    733       // Mask so I can see on the output image where teh footprint is.
    734       for (int i = 0; i < footprint->spans->n; i++) {
    735         pmSpan *sp = footprint->spans->data[i];
    736         for (int j = sp->x0; j <= sp->x1; j++) {
    737           y = sp->y - ym;
    738           x = j - xm;
    739           mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x01;
    740         }
    741       }
    742 
    743       int CRpix_count = 0;     
    744       do {
    745         CRpix_count = 0;
    746         // Zero out my temp images.
    747         for (x = 0; x < binning * dx; x++) {
    748           for (y = 0; y < binning * dy; y++) {
    749             psImageSet(binned,x,y,0.0);
    750             psImageSet(conved,x,y,0.0);
    751             psImageSet(edges,x/binning,y/ binning,0.0);
    752           }
    753         }     
    754         // Make subsampled image. Maybe this should be called "unbinned" or something
    755         for (x = 0; x < binning * dx; x++) {
    756           for (y = 0; y < binning * dy; y++) {
    757             psImageSet(binned,x,y,psImageGet(mypix,x / binning,y / binning));
    758           }
    759         }
    760         // Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero
    761         for (x = 1; x < dx - 1; x++) {
    762           for (y = 1; y < dy - 1; y++) {
    763             psImageSet(conved,x,y,psImageGet(binned,x,y) - 0.25 *
    764                        (psImageGet(binned,x-1,y) + psImageGet(binned,x+1,y) +
    765                         psImageGet(binned,x,y-1) + psImageGet(binned,x,y+1)));
    766             if (psImageGet(conved,x,y) < 0.0) {
    767               psImageSet(conved,x,y,0.0);
    768             }
    769           }
    770         }
    771         // Create an edge map by rebinning
    772         for (x = 0; x < binning * dx; x++) {
    773           for (y = 0; y < binning * dy; y++) {
    774             psImageSet(edges,x / binning, y / binning,
    775                        psImageGet(edges, x / binning, y / binning) +
    776                        psImageGet(conved,x,y));
    777           }
    778         }
    779         // Modify my mask if we're above the significance threshold
    780         for (x = 0; x < dx; x++) {
    781           for (y = 0; y < dy; y++) {
    782             if ( psImageGet(edges,x,y) / (binning * sqrt(noise_factor * psImageGet(myvar,x,y))) > sigma_thresh ) {
    783               if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x01) {
    784                 mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x40;
    785                 CRpix_count++;
    786               }
    787             }
    788           }
    789         }
    790 
    791         // "Repair" Masked pixels for the next round.
    792         for (x = 1; x < dx - 1; x++) {
    793           for (y = 1; y < dy - 1; y++) {
    794             if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
    795               psImageSet(mypix,x,y,0.25 *
    796                          (psImageGet(mypix,x-1,y) + psImageGet(mypix,x+1,y) +
    797                           psImageGet(mypix,x,y-1) + psImageGet(mypix,x,y+1)));
    798             }
    799           }
    800         }
    801        
    802 /*      if ((psTraceGetLevel("psphot.czw") >= 2)&&(abs(xm - 2770) < 5)&&(abs(ym - 2581) < 5)&&(iteration == 0)) { */
    803 /*        psTrace("psphot.czw",2,"TRACEMOTRON %d %d %d %d %d\n",xm,ym,dx,dy,iteration); */
    804 /*        psphotSaveImage (NULL, mypix,   "czw.pix.fits"); */
    805 /*        psphotSaveImage (NULL, myvar,   "czw.var.fits"); */
    806 /*        psphotSaveImage (NULL, binned,  "czw.binn.fits"); */
    807 /*        psphotSaveImage (NULL, conved,  "czw.conv.fits"); */
    808 /*        psphotSaveImage (NULL, edges,   "czw.edge.fits"); */
    809 /*        psphotSaveImage (NULL, mymask,  "czw.mask.fits"); */
    810          
    811 /*      } */
    812        
    813         psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration,CRpix_count);
    814         iteration++;
    815       } while ((iteration < max_iter)&&(CRpix_count > 0));
    816 
    817       // A solitary masked pixel is likely a lie. Remove those.
    818       for (x = 0; x < dx; x++) {
    819         for (y = 0; y < dy; y++) {
    820           if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
    821             if ((x-1 >= 0)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) {
    822               continue;
    823             }
    824             if ((y-1 >= 0)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) {
    825               continue;
    826             }
    827             if ((x+1 < dx)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) {
    828               continue;
    829             }
    830             if ((y+1 < dy)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) {
    831               continue;
    832             }
    833             mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40;
    834           }
    835         }
    836       }
    837              
    838 
    839      
    840       for (x = 0; x < dx; x++) {
    841         for (y = 0; y < dy; y++) {
    842           if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
    843             mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ym+mask->row0][x+xm+mask->col0] |= maskVal;
    844           }
    845 /*        // Look at the footprint */
    846 /*        if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x01) { */
    847 /*          mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ym+mask->row0][x+xm+mask->col0] |= 0x04; */
    848 /*        } */
    849 /*        // Look at the bounding box. */
    850 /*        mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ym+mask->row0][x+xm+mask->col0] |= 0x08; */
    851         }
    852       }
    853      
    854 
    855       psFree(mypix);
    856       psFree(myvar);
    857       psFree(binned);
    858       psFree(conved);
    859       psFree(edges);
    860       psFree(mymask);
    861      
    862     }
    863     else {
    864       psTrace("psphot.czw", 2, "BBox does not match npix: %d %d",(xM - xm) * (yM - ym),footprint->npix);
    865     }   
    866   }
    867   else {
    868     psTrace("psphot.czw", 2, "No valid footprint for source (%g %g)",peak->xf,peak->yf);
    869   }
    870   return(true);
    871 }
    872 
    873      
    874        
Note: See TracChangeset for help on using the changeset viewer.