IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27695


Ignore:
Timestamp:
Apr 15, 2010, 12:45:18 PM (16 years ago)
Author:
Paul Price
Message:

Catch case n=0 to prevent errors getting raised which make psphot fail.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotSourceSize.c

    r27657 r27695  
    5050    // loop over the available readouts
    5151    for (int i = 0; i < num; i++) {
    52         if (i == chisqNum) continue; // skip chisq image
    53         if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i, recipe, getPSFsize)) {
     52        if (i == chisqNum) continue; // skip chisq image
     53        if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i, recipe, getPSFsize)) {
    5454            psError (PSPHOT_ERR_CONFIG, false, "failed on source size analysis for PSPHOT.INPUT entry %d", i);
    55             return false;
    56         }
     55            return false;
     56        }
    5757    }
    5858    return true;
    … …  
    8181
    8282    if (!sources->n) {
    83         psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
    84         return true;
     83        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
     84        return true;
    8585    }
    8686
    … …  
    140140    // XXX this should only be done on the first pass (ie, if we have newSources or allSources?)
    141141    if (getPSFsize) {
    142         psphotSourceSizePSF (&options, sources, psf);
     142        psphotSourceSizePSF (&options, sources, psf);
    143143    }
    144144
    … …  
    173173    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
    174174
     175    int num = 0;                        // Number of sources measured
    175176    for (int i = 0; i < sources->n; i++) {
    176177        pmSource *source = sources->data[i];
    177178        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
     179        num++;
    178180
    179181        // replace object in image
    … …  
    200202        psVectorAppend (Ap, dMag);
    201203        psVectorAppend (ApErr, source->errMag);
     204    }
     205    if (num == 0) {
     206        // Not raising an error, because errors aren't being checked elsewhere in this function
     207        psFree(Ap);
     208        psFree(ApErr);
     209        return false;
    202210    }
    203211
    … …  
    258266            continue;
    259267        }
    260         // psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
     268        // psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
    261269    }
    262270
    … …  
    313321        }
    314322
    315         // convert to Mmaj, Mmin:
     323        // convert to Mmaj, Mmin:
    316324        psF32 Mxx = source->moments->Mxx;
    317325        psF32 Myy = source->moments->Myy;
    … …  
    338346        float dMag = source->psfMag - apMag;
    339347
    340         // set nSigma to include both systematic and poisson error terms
    341         // XXX the 'poisson error' contribution for size is probably wrong...
     348        // set nSigma to include both systematic and poisson error terms
     349        // XXX the 'poisson error' contribution for size is probably wrong...
    342350        float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);
    343         float nSigmaMXX = (Mxx - psfClump->X) / hypot(psfClump->dX, psfClump->X*psfClump->X*source->errMag);
    344         float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag);
    345 
    346         // partially-masked sources are more likely to be mis-measured PSFs
    347         float sizeBias = 1.0;
    348         if (source->pixWeight < 0.9) {
    349             sizeBias = 3.0;
    350         }
    351 
    352         float minMxx = psfClump->X - sizeBias*options->nSigmaMoments*psfClump->dX;
    353         float minMyy = psfClump->Y - sizeBias*options->nSigmaMoments*psfClump->dY;
    354 
    355         // include MAG, MXX, and MYY?
     351        float nSigmaMXX = (Mxx - psfClump->X) / hypot(psfClump->dX, psfClump->X*psfClump->X*source->errMag);
     352        float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag);
     353
     354        // partially-masked sources are more likely to be mis-measured PSFs
     355        float sizeBias = 1.0;
     356        if (source->pixWeight < 0.9) {
     357            sizeBias = 3.0;
     358        }
     359
     360        float minMxx = psfClump->X - sizeBias*options->nSigmaMoments*psfClump->dX;
     361        float minMyy = psfClump->Y - sizeBias*options->nSigmaMoments*psfClump->dY;
     362
     363        // include MAG, MXX, and MYY?
    356364        source->extNsigma = nSigmaMAG;
    357365
    358         // notes to clarify the source size classification rules:
    359         // * a defect should be functionally equivalent to a cosmic ray
    360         // * CR & defect should have a faintess limit (min S/N)
    361         // * SAT stars should not be faint, but defects may?
     366        // notes to clarify the source size classification rules:
     367        // * a defect should be functionally equivalent to a cosmic ray
     368        // * CR & defect should have a faintess limit (min S/N)
     369        // * SAT stars should not be faint, but defects may?
    362370
    363371        // Anything within this region is a probably PSF-like object. Saturated stars may land
    … …  
    365373        bool isPSF = (fabs(nSigmaMAG) < options->nSigmaApResid) && (fabs(nSigmaMXX) < sizeBias*options->nSigmaMoments) && (fabs(nSigmaMYY) < sizeBias*options->nSigmaMoments);
    366374        if (isPSF) {
    367           psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
    368                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
    369                   options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    370           source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    371           Npsf ++;
    372           continue;
     375          psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
     376                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
     377                  options->nSigmaApResid,sizeBias*options->nSigmaMoments);
     378          source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     379          Npsf ++;
     380          continue;
    373381        }
    374382
    … …  
    376384        // Defects may also be marked as SATSTAR -- XXX deactivate this flag?
    377385        // XXX this rule is not great
    378         // XXX only accept brightish detections as CRs
    379         // (nSigmaMAG < -options->nSigmaApResid) ||
    380         bool isCR = isCR = (source->errMag < 1.0 / SIZE_SN_LIM) && ((Mxx < minMxx) || (Myy < minMyy));
    381         if (isCR) {
    382             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",
    383                     source->peak->xf,source->peak->yf,source->pixWeight,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
    384                     options->nSigmaApResid,sizeBias*options->nSigmaMoments);
     386        // XXX only accept brightish detections as CRs
     387        // (nSigmaMAG < -options->nSigmaApResid) ||
     388        bool isCR = isCR = (source->errMag < 1.0 / SIZE_SN_LIM) && ((Mxx < minMxx) || (Myy < minMyy));
     389        if (isCR) {
     390            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",
     391                    source->peak->xf,source->peak->yf,source->pixWeight,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
     392                    options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    385393            source->mode |= PM_SOURCE_MODE_DEFECT;
    386             source->tmpFlags |= PM_SOURCE_TMPF_SIZE_CR_CANDIDATE;
     394            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_CR_CANDIDATE;
    387395            Ncr ++;
    388396            continue;
    … …  
    392400        // just large saturated regions.
    393401        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    394             psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g SAT\t%g %g\n",
    395                     source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
    396                     options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    397             source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     402            psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g SAT\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);
     405            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    398406            Nsat ++;
    399407            continue;
    … …  
    403411        bool isEXT = (nSigmaMAG > options->nSigmaApResid) || (Mxx > minMxx) || (Myy > minMyy);
    404412        if (isEXT) {
    405           psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Ext\t%g %g\n",
    406                   source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
    407                   options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    408          
     413          psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Ext\t%g %g\n",
     414                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
     415                  options->nSigmaApResid,sizeBias*options->nSigmaMoments);
     416
    409417            source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
    410             source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     418            source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    411419            Next ++;
    412420            continue;
    413421        }
    414         psTrace("psphotSourceClassRegion.MISS",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Unk\t%g %g\n",
    415                 source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
    416                 options->nSigmaApResid,sizeBias*options->nSigmaMoments);
    417        
    418         // sources that reach here are probably too faint for a reasonable source size measurement
     422        psTrace("psphotSourceClassRegion.MISS",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Unk\t%g %g\n",
     423                source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG,
     424                options->nSigmaApResid,sizeBias*options->nSigmaMoments);
     425
     426        // sources that reach here are probably too faint for a reasonable source size measurement
    419427        // psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigmaMAG);
    420         source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     428        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    421429        Nmiss ++;
    422430    }
    … …  
    450458
    451459        // skip unless this source is thought to be a cosmic ray.  flag the detection and mask the pixels
    452         // XXX this may be degenerate with the above test
    453         if (!(source->mode & PM_SOURCE_MODE_DEFECT)) continue;
     460        // XXX this may be degenerate with the above test
     461        if (!(source->mode & PM_SOURCE_MODE_DEFECT)) continue;
    454462
    455463        // Integer position of peak
    … …  
    463471            continue;
    464472        }
    465        
    466         // XXX for testing, only CRMASK a single source:
    467         if (options->xtest && (fabs(source->peak->xf - options->xtest) > 5)) continue;
    468         if (options->ytest && (fabs(source->peak->yf - options->ytest) > 5)) continue;
     473
     474        // XXX for testing, only CRMASK a single source:
     475        if (options->xtest && (fabs(source->peak->xf - options->xtest) > 5)) continue;
     476        if (options->ytest && (fabs(source->peak->yf - options->ytest) > 5)) continue;
    469477
    470478        // replace object in image
    … …  
    473481        }
    474482
    475         // XXX this is running slowly and is too agressive, but it more-or-less works
    476         psTrace("psphot", 6, "mask cosmic ray at %f, %f\n", source->peak->xf, source->peak->yf);
    477         if (options->apply) {
    478             psphotMaskCosmicRay(readout, source, options->crMask);
    479         } else {
    480             source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    481         }
    482         nMasked ++;
     483        // XXX this is running slowly and is too agressive, but it more-or-less works
     484        psTrace("psphot", 6, "mask cosmic ray at %f, %f\n", source->peak->xf, source->peak->yf);
     485        if (options->apply) {
     486            psphotMaskCosmicRay(readout, source, options->crMask);
     487        } else {
     488            source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     489        }
     490        nMasked ++;
    483491
    484492        // re-subtract the object, leave local sky
    485493        pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
    486494    }
    487    
     495
    488496    // now that we have masked pixels associated with CRs, we can grow the mask
    489497    if (options->grow > 0) {
    … …  
    503511    // XXX test : save the mask image
    504512    if (0) {
    505         psphotSaveImage (NULL, readout->mask,   "mask.fits");
     513        psphotSaveImage (NULL, readout->mask,   "mask.fits");
    506514    }
    507515
    … …  
    540548    psImage *image= readout->image;
    541549    psImage *variance = readout->variance;
    542      
     550
    543551    int binning = 2;
    544552    float sigma_thresh = 3.0;
    … …  
    553561    psImage *edges  = psImageAlloc(dx,dy,image->type.type);
    554562    psImage *mymask = psImageAlloc(dx,dy,PS_TYPE_IMAGE_MASK);
    555      
     563
    556564    // Load my copy of things.
    557565    for (int y = 0; y < dy; y++) {
    558         for (int x = 0; x < dx; x++) {
    559             mypix->data.F32[y][x] = image->data.F32[y+ys][x+xs];
    560             myvar->data.F32[y][x] = variance->data.F32[y+ys][x+xs];
    561             mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0x00;
    562         }
     566        for (int x = 0; x < dx; x++) {
     567            mypix->data.F32[y][x] = image->data.F32[y+ys][x+xs];
     568            myvar->data.F32[y][x] = variance->data.F32[y+ys][x+xs];
     569            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0x00;
     570        }
    563571    }
    564572    // Mask so I can see on the output image where the footprint is.
    565573    for (int i = 0; i < footprint->spans->n; i++) {
    566         pmSpan *sp = footprint->spans->data[i];
    567         for (int j = sp->x0; j <= sp->x1; j++) {
    568             int y = sp->y - ys;
    569             int x = j - xs;
    570             mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x01;
    571         }
     574        pmSpan *sp = footprint->spans->data[i];
     575        for (int j = sp->x0; j <= sp->x1; j++) {
     576            int y = sp->y - ys;
     577            int x = j - xs;
     578            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x01;
     579        }
    572580    }
    573581
    574582    int nCRpix = 1; // force at least one pass...
    575583    for (int iteration = 0; (iteration < max_iter) && (nCRpix > 0); iteration++) {
    576         nCRpix = 0;
    577         psImageInit (binned, 0.0);
    578         psImageInit (conved, 0.0);
    579         psImageInit (edges, 0.0);
    580 
    581         // Make subsampled image. Maybe this should be called "unbinned" or something
    582         for (int y = 0; y < binning * dy; y++) {
    583             int yraw = y / binning;
    584             for (int x = 0; x < binning * dx; x++) {
    585                 int xraw = x / binning;
    586                 binned->data.F32[y][x] = mypix->data.F32[yraw][xraw];
    587             }
    588         }
    589 
    590         // Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero
    591         for (int y = 1; y < binning * dy - 1; y++) {
    592             for (int x = 1; x < binning * dx - 1; x++) {
    593                 float value = binned->data.F32[y][x] - 0.25 *
    594                     (binned->data.F32[y+0][x-1] + binned->data.F32[y+0][x+1] +
    595                      binned->data.F32[y-1][x+0] + binned->data.F32[y+1][x+0]);
    596                 value = PS_MAX(0.0, value);
    597 
    598                 conved->data.F32[y][x] = value;
    599             }
    600         }
    601 
    602         // Create an edge map by rebinning
    603         for (int y = 0; y < binning * dy; y++) {
    604             int yraw = y / binning;
    605             for (int x = 0; x < binning * dx; x++) {
    606                 int xraw = x / binning;
    607                 edges->data.F32[yraw][xraw] += conved->data.F32[y][x];
    608             }
    609         }
    610 
    611         // coordinate of peak in subimage pixels:
    612         int xPeak = peak->x - xs;
    613         int yPeak = peak->y - ys;
    614 
    615         // Modify my mask if we're above the significance threshold, but only for connected pixels
    616         nCRpix = psphotMaskCosmicRayConnected (xPeak, yPeak, mymask, myvar, edges, binning, sigma_thresh);
    617        
     584        nCRpix = 0;
     585        psImageInit (binned, 0.0);
     586        psImageInit (conved, 0.0);
     587        psImageInit (edges, 0.0);
     588
     589        // Make subsampled image. Maybe this should be called "unbinned" or something
     590        for (int y = 0; y < binning * dy; y++) {
     591            int yraw = y / binning;
     592            for (int x = 0; x < binning * dx; x++) {
     593                int xraw = x / binning;
     594                binned->data.F32[y][x] = mypix->data.F32[yraw][xraw];
     595            }
     596        }
     597
     598        // Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero
     599        for (int y = 1; y < binning * dy - 1; y++) {
     600            for (int x = 1; x < binning * dx - 1; x++) {
     601                float value = binned->data.F32[y][x] - 0.25 *
     602                    (binned->data.F32[y+0][x-1] + binned->data.F32[y+0][x+1] +
     603                     binned->data.F32[y-1][x+0] + binned->data.F32[y+1][x+0]);
     604                value = PS_MAX(0.0, value);
     605
     606                conved->data.F32[y][x] = value;
     607            }
     608        }
     609
     610        // Create an edge map by rebinning
     611        for (int y = 0; y < binning * dy; y++) {
     612            int yraw = y / binning;
     613            for (int x = 0; x < binning * dx; x++) {
     614                int xraw = x / binning;
     615                edges->data.F32[yraw][xraw] += conved->data.F32[y][x];
     616            }
     617        }
     618
     619        // coordinate of peak in subimage pixels:
     620        int xPeak = peak->x - xs;
     621        int yPeak = peak->y - ys;
     622
     623        // Modify my mask if we're above the significance threshold, but only for connected pixels
     624        nCRpix = psphotMaskCosmicRayConnected (xPeak, yPeak, mymask, myvar, edges, binning, sigma_thresh);
     625
    618626# if DUMPPICS
    619         psphotSaveImage (NULL, mypix,   "crmask.pix.fits");
     627        psphotSaveImage (NULL, mypix,   "crmask.pix.fits");
    620628# endif
    621        
     629
    622630// XXX do not repair the pixels in isophot version
    623631# if 0
    624         // "Repair" Masked pixels for the next round.
    625         for (int y = 1; y < dy - 1; y++) {
    626             for (int x = 1; x < dx - 1; x++) {
    627                 if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) {
    628                     myfix->data.F32[y][x] = mypix->data.F32[y][x];
    629                     continue;
    630                 }
    631                 myfix->data.F32[y][x] = 0.25 *
    632                     (mypix->data.F32[y+0][x-1] + mypix->data.F32[y+0][x+1] +
    633                      mypix->data.F32[y-1][x+0] + mypix->data.F32[y+1][x+0]);
    634             }
    635         }
    636        
    637         // "Repair" Masked pixels for the next round.
    638         for (int y = 1; y < dy - 1; y++) {
    639             for (int x = 1; x < dx - 1; x++) {
    640                 mypix->data.F32[y][x] = myfix->data.F32[y][x];
    641             }
    642         }
     632        // "Repair" Masked pixels for the next round.
     633        for (int y = 1; y < dy - 1; y++) {
     634            for (int x = 1; x < dx - 1; x++) {
     635                if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) {
     636                    myfix->data.F32[y][x] = mypix->data.F32[y][x];
     637                    continue;
     638                }
     639                myfix->data.F32[y][x] = 0.25 *
     640                    (mypix->data.F32[y+0][x-1] + mypix->data.F32[y+0][x+1] +
     641                     mypix->data.F32[y-1][x+0] + mypix->data.F32[y+1][x+0]);
     642            }
     643        }
     644
     645        // "Repair" Masked pixels for the next round.
     646        for (int y = 1; y < dy - 1; y++) {
     647            for (int x = 1; x < dx - 1; x++) {
     648                mypix->data.F32[y][x] = myfix->data.F32[y][x];
     649            }
     650        }
    643651# endif
    644652
    645653# if DUMPPICS
    646         fprintf (stderr, "CRMASK %d %d %d %d %d\n", xs, ys, dx, dy, iteration);
    647         psphotSaveImage (NULL, mypix,   "crmask.fix.fits");
    648         psphotSaveImage (NULL, myvar,   "crmask.var.fits");
    649         psphotSaveImage (NULL, binned,  "crmask.binn.fits");
    650         psphotSaveImage (NULL, conved,  "crmask.conv.fits");
    651         psphotSaveImage (NULL, edges,   "crmask.edge.fits");
    652         psphotSaveImage (NULL, mymask,  "crmask.mask.fits");
     654        fprintf (stderr, "CRMASK %d %d %d %d %d\n", xs, ys, dx, dy, iteration);
     655        psphotSaveImage (NULL, mypix,   "crmask.fix.fits");
     656        psphotSaveImage (NULL, myvar,   "crmask.var.fits");
     657        psphotSaveImage (NULL, binned,  "crmask.binn.fits");
     658        psphotSaveImage (NULL, conved,  "crmask.conv.fits");
     659        psphotSaveImage (NULL, edges,   "crmask.edge.fits");
     660        psphotSaveImage (NULL, mymask,  "crmask.mask.fits");
    653661# endif
    654         psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration, nCRpix);
     662        psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration, nCRpix);
    655663    }
    656664
    … …  
    659667    // XXX can't we use nCRpix == 1 to test for these?
    660668    for (int x = 0; x < dx; x++) {
    661         for (int y = 0; y < dy; y++) {
    662             if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) continue;
    663             if ((x-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) {
    664                 continue;
    665             }
    666             if ((y-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) {
    667                 continue;
    668             }
    669             if ((x+1 < dx) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) {
    670                 continue;
    671             }
    672             if ((y+1 < dy) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) {
    673                 continue;
    674             }
    675             mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40;
    676         }
     669        for (int y = 0; y < dy; y++) {
     670            if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) continue;
     671            if ((x-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) {
     672                continue;
     673            }
     674            if ((y-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) {
     675                continue;
     676            }
     677            if ((x+1 < dx) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) {
     678                continue;
     679            }
     680            if ((y+1 < dy) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) {
     681                continue;
     682            }
     683            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40;
     684        }
    677685    }
    678686# endif
    … …  
    681689    nCRpix = 0;
    682690    for (int x = 0; x < dx; x++) {
    683         for (int y = 0; y < dy; y++) {
    684             if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
    685                 mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ys+mask->row0][x+xs+mask->col0] |= maskVal;
    686                 nCRpix ++;
    687             }
    688         }
    689     }
    690      
     691        for (int y = 0; y < dy; y++) {
     692            if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
     693                mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ys+mask->row0][x+xs+mask->col0] |= maskVal;
     694                nCRpix ++;
     695            }
     696        }
     697    }
     698
    691699    // XXX if we decide this REALLY is a cosmic ray, set the CR_LIMIT bit
    692700    if (nCRpix > 1) {
    693         source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    694         source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     701        source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     702        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    695703    }
    696704    // fprintf (stderr, "CRMASK %d %d %d %d %d\n", peak->x, peak->y, dx, dy, nCRpix);
    … …  
    703711    psFree(edges);
    704712    psFree(mymask);
    705      
     713
    706714    return true;
    707715}
    … …  
    711719    for (int i = 0; i < sources->n; i++) {
    712720        pmSource *source = sources->data[i];
    713         pmPeak *peak = source->peak;
    714         pmFootprint *footprint = peak->footprint;
    715         if (!footprint) continue;
    716         for (int j = 0; j < footprint->spans->n; j++) {
    717             pmSpan *sp = footprint->spans->data[j];
    718             psAssert (sp, "missing span");
    719         }
     721        pmPeak *peak = source->peak;
     722        pmFootprint *footprint = peak->footprint;
     723        if (!footprint) continue;
     724        for (int j = 0; j < footprint->spans->n; j++) {
     725            pmSpan *sp = footprint->spans->data[j];
     726            psAssert (sp, "missing span");
     727        }
    720728    }
    721729    return true;
    … …  
    741749    if (!footprint) {
    742750      psTrace("psphot.czw",2,"Using isophot CR mask code.");
    743      
     751
    744752        // if we have not footprint, use the old code to mask by isophot
    745753        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    … …  
    749757    if (!footprint->spans) {
    750758      psTrace("psphot.czw",2,"Using isophot CR mask code.");
    751      
     759
    752760        // if we have no footprint, use the old code to mask by isophot
    753761        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    … …  
    756764    psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    757765    // mask all of the pixels covered by the spans of the footprint
    758     for (int j = 1; j < footprint->spans->n; j++) { 
    759         pmSpan *span1 = footprint->spans->data[j];
    760 
    761         int iy = span1->y;
    762         int xs = span1->x0;
    763         int xe = span1->x1;
    764 
    765         for (int ix = xs; ix < xe; ix++) {
    766             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    767         }
    768     } 
     766    for (int j = 1; j < footprint->spans->n; j++) {
     767        pmSpan *span1 = footprint->spans->data[j];
     768
     769        int iy = span1->y;
     770        int xs = span1->x0;
     771        int xe = span1->x1;
     772
     773        for (int ix = xs; ix < xe; ix++) {
     774            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     775        }
     776    }
    769777    return true;
    770778}
    … …  
    772780# define VERBOSE 0
    773781int psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh) {
    774    
     782
    775783    int xLo, xRo;
    776784    int nCRpix = 0;
    777785
    778786    float noise_factor = 5.0 / 4.0;  // Intrinsic to the Laplacian making noise spikes spikier.
    779    
     787
    780788    // mark the pixels in this row to the left, then the right. stay within footprint
    781789    int xL = xPeak; // find the range of valid pixels in this row
    782790    int xR = xPeak;
    783791    for (int ix = xPeak; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] & 0x01); ix--) {
    784         float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]);
    785         float value = edges->data.F32[yPeak][ix] / noise;
    786         if (value < sigma_thresh ) break;
    787         mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40;
    788         xL = ix;
    789         nCRpix ++;
    790         if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR);
     792        float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]);
     793        float value = edges->data.F32[yPeak][ix] / noise;
     794        if (value < sigma_thresh ) break;
     795        mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40;
     796        xL = ix;
     797        nCRpix ++;
     798        if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR);
    791799    }
    792800    for (int ix = xPeak; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] & 0x01); ix++) {
    793         float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]);
    794         float value = edges->data.F32[yPeak][ix] / noise;
    795         if (value < sigma_thresh ) break;
    796         mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40;
    797         xR = ix;
    798         nCRpix ++;
    799         if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR);
     801        float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]);
     802        float value = edges->data.F32[yPeak][ix] / noise;
     803        if (value < sigma_thresh ) break;
     804        mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40;
     805        xR = ix;
     806        nCRpix ++;
     807        if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR);
    800808    }
    801809    // xL and xR mark the first and last valid pixel in the row
    802810
    803     // for each of the neighboring rows, mark the high pixels if they touch the range xL to xR 
     811    // for each of the neighboring rows, mark the high pixels if they touch the range xL to xR
    804812    xLo = PS_MAX(xL - 1, 0);
    805813    xRo = PS_MIN(xR + 1, mymask->numCols);
    … …  
    808816    for (int iy = yPeak - 1; iy >= 0; iy--) {
    809817
    810         int xLn = -1;
    811         int xRn = -1;
    812         int newPix = 0;
    813 
    814         // mark the pixels in the good range
    815         for (int ix = xLo; ix < xRo; ix++) {
    816             if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint
    817             float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
    818             float value = edges->data.F32[iy][ix] / noise;
    819             if (value < sigma_thresh ) continue;
    820             mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
    821             if (xLn == -1) xLn = ix; // first valid pixel in this row
    822             xRn = ix;                // last valid pixel in this row
    823             nCRpix ++;
    824             newPix ++;
    825             if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
    826         }
    827 
    828         // mark the pixels to the left of the good range
    829         for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) {
    830             float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
    831             float value = edges->data.F32[iy][ix] / noise;
    832             if (value < sigma_thresh ) break;
    833             mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
    834             if (xRn == -1) xRn = ix; // last valid pixel in this row
    835             xLn = ix;
    836             nCRpix ++;
    837             newPix ++;
    838             if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
    839         }
    840 
    841         // mark the pixels to the right of the good range
    842         for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) {
    843             float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
    844             float value = edges->data.F32[iy][ix] / noise;
    845             if (value < sigma_thresh ) break;
    846             mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
    847             if (xLn == -1) xLn = ix; // first valid pixel in this row
    848             xRn = ix;
    849             nCRpix ++;
    850             newPix ++;
    851             if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
    852         }
    853         if (newPix == 0) break;
    854         xLo = PS_MAX(xLn - 1, 0);
    855         xRo = PS_MIN(xRn + 1, mymask->numCols);
     818        int xLn = -1;
     819        int xRn = -1;
     820        int newPix = 0;
     821
     822        // mark the pixels in the good range
     823        for (int ix = xLo; ix < xRo; ix++) {
     824            if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint
     825            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
     826            float value = edges->data.F32[iy][ix] / noise;
     827            if (value < sigma_thresh ) continue;
     828            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
     829            if (xLn == -1) xLn = ix; // first valid pixel in this row
     830            xRn = ix;                // last valid pixel in this row
     831            nCRpix ++;
     832            newPix ++;
     833            if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
     834        }
     835
     836        // mark the pixels to the left of the good range
     837        for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) {
     838            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
     839            float value = edges->data.F32[iy][ix] / noise;
     840            if (value < sigma_thresh ) break;
     841            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
     842            if (xRn == -1) xRn = ix; // last valid pixel in this row
     843            xLn = ix;
     844            nCRpix ++;
     845            newPix ++;
     846            if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
     847        }
     848
     849        // mark the pixels to the right of the good range
     850        for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) {
     851            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
     852            float value = edges->data.F32[iy][ix] / noise;
     853            if (value < sigma_thresh ) break;
     854            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
     855            if (xLn == -1) xLn = ix; // first valid pixel in this row
     856            xRn = ix;
     857            nCRpix ++;
     858            newPix ++;
     859            if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
     860        }
     861        if (newPix == 0) break;
     862        xLo = PS_MAX(xLn - 1, 0);
     863        xRo = PS_MIN(xRn + 1, mymask->numCols);
    856864    }
    857865
    … …  
    862870    for (int iy = yPeak + 1; iy < mymask->numRows; iy++) {
    863871
    864         int xLn = -1;
    865         int xRn = -1;
    866         int newPix = 0;
    867 
    868         // mark the pixels in the good range
    869         for (int ix = xLo; ix < xRo; ix++) {
    870             if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint
    871             float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
    872             float value = edges->data.F32[iy][ix] / noise;
    873             if (value < sigma_thresh ) continue;
    874             mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
    875             if (xLn == -1) xLn = ix; // first valid pixel in this row
    876             xRn = ix;                // last valid pixel in this row
    877             nCRpix ++;
    878             newPix ++;
    879             if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
    880         }
    881 
    882         // mark the pixels to the left of the good range
    883         for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) {
    884             float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
    885             float value = edges->data.F32[iy][ix] / noise;
    886             if (value < sigma_thresh ) break;
    887             mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
    888             if (xRn == -1) xRn = ix; // last valid pixel in this row
    889             xLn = ix;
    890             nCRpix ++;
    891             newPix ++;
    892             if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
    893         }
    894 
    895         // mark the pixels to the right of the good range
    896         for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) {
    897             float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
    898             float value = edges->data.F32[iy][ix] / noise;
    899             if (value < sigma_thresh ) break;
    900             mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
    901             if (xLn == -1) xLn = ix; // first valid pixel in this row
    902             xRn = ix;
    903             nCRpix ++;
    904             newPix ++;
    905             if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
    906         }
    907         if (newPix == 0) break;
    908         xLo = PS_MAX(xLn - 1, 0);
    909         xRo = PS_MIN(xRn + 1, mymask->numCols);
     872        int xLn = -1;
     873        int xRn = -1;
     874        int newPix = 0;
     875
     876        // mark the pixels in the good range
     877        for (int ix = xLo; ix < xRo; ix++) {
     878            if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint
     879            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
     880            float value = edges->data.F32[iy][ix] / noise;
     881            if (value < sigma_thresh ) continue;
     882            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
     883            if (xLn == -1) xLn = ix; // first valid pixel in this row
     884            xRn = ix;                // last valid pixel in this row
     885            nCRpix ++;
     886            newPix ++;
     887            if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
     888        }
     889
     890        // mark the pixels to the left of the good range
     891        for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) {
     892            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
     893            float value = edges->data.F32[iy][ix] / noise;
     894            if (value < sigma_thresh ) break;
     895            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
     896            if (xRn == -1) xRn = ix; // last valid pixel in this row
     897            xLn = ix;
     898            nCRpix ++;
     899            newPix ++;
     900            if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
     901        }
     902
     903        // mark the pixels to the right of the good range
     904        for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) {
     905            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
     906            float value = edges->data.F32[iy][ix] / noise;
     907            if (value < sigma_thresh ) break;
     908            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
     909            if (xLn == -1) xLn = ix; // first valid pixel in this row
     910            xRn = ix;
     911            nCRpix ++;
     912            newPix ++;
     913            if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
     914        }
     915        if (newPix == 0) break;
     916        xLo = PS_MAX(xLn - 1, 0);
     917        xRo = PS_MIN(xRn + 1, mymask->numCols);
    910918    }
    911919
Note: See TracChangeset for help on using the changeset viewer.