IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26617


Ignore:
Timestamp:
Jan 15, 2010, 3:21:15 PM (16 years ago)
Author:
watersc1
Message:

I don't know why this wasn't there already.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotSourceSize.c

    r25852 r26617  
    2222bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
    2323
     24bool psphotMaskCosmicRayCZW (pmReadout *readout, pmSource *source, psImageMaskType maskVal,psphotSourceSizeOptions *options);
    2425// we need to call this function after sources have been fitted to the PSF model and
    2526// subtracted.  To determine the CR-nature, this function examines the 9 pixels in the 3x3
    … …  
    7071
    7172    // We are using the value psfMag - 2.5*log10(moment.Sum) as a measure of the extendedness
    72     // of and object.  We need to model this distribution for the PSF stars before we can test
     73    // of an object.  We need to model this distribution for the PSF stars before we can test
    7374    // the significance for a specific object
    7475    // XXX move this to the code that generates the PSF?
    … …  
    103104    pmFootprint *footprint = peak->footprint;
    104105    if (!footprint) {
     106      psTrace("psphot.czw",2,"Using isophot CR mask code.");
     107     
    105108        // if we have not footprint, use the old code to mask by isophot
    106109        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    … …  
    109112
    110113    if (!footprint->spans) {
     114      psTrace("psphot.czw",2,"Using isophot CR mask code.");
     115     
    111116        // if we have no footprint, use the old code to mask by isophot
    112117        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    113118        return true;
    114119    }
    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     }
     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/*     } */
    128133    return true;
    129134}
     135
     136 
     137
    130138
    131139bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
    … …  
    373381        // Anything within this region is a probably PSF-like object. Saturated stars may land
    374382        // in this region, but are detected elsewhere on the basis of their peak value.
    375         bool isPSF = (fabs(nSigma) < options->nSigmaApResid) && (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) && (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY);
     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));
    376386        if (isPSF) {
    377             Npsf ++;
    378             continue;
     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;
    379392        }
    380393
    … …  
    383396        // XXX this rule is not great
    384397        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);
    385402            source->mode |= PM_SOURCE_MODE_DEFECT;
    386403            Ncr ++;
    … …  
    391408        // just large saturated regions.
    392409        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         
    393414            Nsat ++;
    394415            continue;
    … …  
    398419        bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y));
    399420        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         
    400425            source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
    401426            Next ++;
    402427            continue;
    403428        }
    404 
     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       
    405433        psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma);
    406434        Nmiss ++;
    … …  
    493521        pmSource *source = sources->data[i];
    494522
    495         // skip source if it was already measured
    496         if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
    497             psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
    498             continue;
    499         }
    500 
    501         // source must have been subtracted
    502         if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
    503             source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
    504             psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
    505             continue;
    506         }
    507 
    508         psF32 **resid  = source->pixels->data.F32;
    509         psF32 **variance = source->variance->data.F32;
    510         psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     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;
    511541
    512542        // Integer position of peak
    … …  
    518548            yPeak < 1 || yPeak > source->pixels->numRows - 2) {
    519549            psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
    520             continue;
    521         }
    522 
    523         // Skip sources with masked pixels.  These are mostly caught as DEFECT
    524         bool keep = true;
    525         for (int iy = -1; (iy <= +1) && keep; iy++) {
    526             for (int ix = -1; (ix <= +1) && keep; ix++) {
    527                 if (mask[yPeak+iy][xPeak+ix] & options->maskVal) {
    528                     keep = false;
    529                 }
    530             }
    531         }
    532         if (!keep) {
    533             psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n");
    534             continue;
    535         }
    536 
    537         // Compare the central pixel with those on either side, for the four possible lines through it.
    538 
    539         // Soften variances (add systematic error)
    540         float softening = options->soft * PS_SQR(source->peak->flux); // Softening for variances
    541 
    542         // Across the middle: y = 0
    543         float cX = 2*resid[yPeak][xPeak]   - resid[yPeak+0][xPeak-1]  - resid[yPeak+0][xPeak+1];
    544         float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1];
    545         float nX = cX / sqrtf(dcX + softening);
    546 
    547         // Up the centre: x = 0
    548         float cY = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak+0]  - resid[yPeak+1][xPeak+0];
    549         float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0];
    550         float nY = cY / sqrtf(dcY + softening);
    551 
    552         // Diagonal: x = y
    553         float cL = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak-1]  - resid[yPeak+1][xPeak+1];
    554         float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1];
    555         float nL = cL / sqrtf(dcL + softening);
    556 
    557         // Diagonal: x = - y
    558         float cR = 2*resid[yPeak][xPeak]   - resid[yPeak+1][xPeak-1]  - resid[yPeak-1][xPeak+1];
    559         float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1];
    560         float nR = cR / sqrtf(dcR + softening);
    561 
    562         // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
    563         // Ndof = 4 ? (four measurements, no free parameters)
    564         // XXX this value is going to be biased low because of systematic errors.
    565         // we need to calibrate it somehow
    566         // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq);
    567 
    568         // not strictly accurate: overcounts the chisq contribution from the center pixel (by
    569         // factor of 4); also biases a bit low if any pixels are masked
    570         // XXX I am not sure I want to keep this value...
    571         source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR);
    572 
    573         float fCR = 0.0;
    574         int nCR = 0;
    575         if (nX > 0.0) {
    576             fCR += nX;
    577             nCR ++;
    578         }
    579         if (nY > 0.0) {
    580             fCR += nY;
    581             nCR ++;
    582         }
    583         if (nL > 0.0) {
    584             fCR += nL;
    585             nCR ++;
    586         }
    587         if (nR > 0.0) {
    588             fCR += nR;
    589             nCR ++;
    590         }
    591         source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
    592         source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    593 
    594         if (!isfinite(source->crNsigma)) {
    595             continue;
    596         }
    597 
     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       
    598633        // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
    599         if (source->crNsigma > options->nSigmaCR) {
    600             source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    601             // XXX still testing... : psphotMaskCosmicRay (readout->mask, source, maskVal, crMask);
    602             // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask);
    603         }
    604     }
    605 
     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   
    606657    // now that we have masked pixels associated with CRs, we can grow the mask
    607658    if (options->grow > 0) {
    … …  
    616667        readout->mask = newMask;
    617668    }
     669
     670    if (psTraceGetLevel("psphot.czw") >= 2) {
     671      psphotSaveImage (NULL, readout->mask,   "mask.fits");
     672    }
    618673    return true;
    619674}
     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
     681bool 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.