Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 20878) +++ trunk/psphot/src/psphot.h (revision 20938) @@ -165,5 +165,5 @@ bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources); bool psphotVisualShowFlags (psArray *sources); -bool psphotVisualShowSourceSize (psArray *sources); +bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources); bool psphotVisualShowResidualImage (pmReadout *readout); bool psphotVisualPlotApResid (psArray *sources); Index: trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinear.c (revision 20878) +++ trunk/psphot/src/psphotFitSourcesLinear.c (revision 20938) @@ -68,5 +68,5 @@ // if (source->type == PM_SOURCE_TYPE_STAR && - // source->mode & PM_SOURCE_MODE_SATSTAR) continue; + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; if (final) { Index: trunk/psphot/src/psphotSourceSize.c =================================================================== --- trunk/psphot/src/psphotSourceSize.c (revision 20878) +++ trunk/psphot/src/psphotSourceSize.c (revision 20938) @@ -4,4 +4,7 @@ static float psphotModelContour(const psImage *image, const psImage *weight, const psImage *mask, psMaskType maskVal, const pmModel *model, float Ro); + +bool psphotMaskCosmicRay_Old (pmSource *source, psMaskType maskVal, psMaskType crMask); +bool psphotMaskCosmicRay_New (psImage *mask, pmSource *source, psMaskType maskVal, psMaskType crMask); // we need to call this function after sources have been fitted to the PSF model and @@ -157,69 +160,6 @@ // this source is thought to be a cosmic ray. flag the detection and mask the pixels if (source->crNsigma > CR_NSIGMA_LIMIT) { - source->mode |= PM_SOURCE_MODE_CR_LIMIT; - pmPeak *peak = source->peak; - psAssert (peak, "NULL peak"); - - // replace the source flux - pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - source->mode &= ~PM_SOURCE_MODE_SUBTRACTED; - - psImage *mask = source->maskView; - psImage *pixels = source->pixels; - psImage *weight = source->weight; - - // XXX This should be a recipe variable - # define SN_LIMIT 5.0 - - int xo = peak->x - pixels->col0; - int yo = peak->y - pixels->row0; - - // mark the pixels in this row to the left, then the right - for (int ix = xo; ix >= 0; ix--) { - float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]); - if (SN > SN_LIMIT) { - mask->data.U8[yo][ix] |= crMask; - } - } - for (int ix = xo + 1; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]); - if (SN > SN_LIMIT) { - mask->data.U8[yo][ix] |= crMask; - } - } - - // for each of the neighboring rows, mark the high pixels if they have a marked neighbor - // first go up: - for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) { - // mark the pixels in this row to the left, then the right - for (int ix = 0; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]); - if (SN < SN_LIMIT) continue; - - bool valid = false; - valid |= (mask->data.U8[iy+1][ix] & crMask); - valid |= (ix > 0) ? (mask->data.U8[iy+1][ix-1] & crMask) : 0; - valid |= (ix <= mask->numCols) ? (mask->data.U8[iy+1][ix+1] & crMask) : 0; - - if (!valid) continue; - mask->data.U8[iy][ix] |= crMask; - } - } - // next go down: - for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) { - // mark the pixels in this row to the left, then the right - for (int ix = 0; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]); - if (SN < SN_LIMIT) continue; - - bool valid = false; - valid |= (mask->data.U8[iy-1][ix] & crMask); - valid |= (ix > 0) ? (mask->data.U8[iy-1][ix-1] & crMask) : 0; - valid |= (ix <= mask->numCols) ? (mask->data.U8[iy-1][ix+1] & crMask) : 0; - - if (!valid) continue; - mask->data.U8[iy][ix] |= crMask; - } - } + // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask); + psphotMaskCosmicRay_Old (source, maskVal, crMask); } } @@ -242,5 +182,5 @@ psphotVisualPlotSourceSize (sources); - psphotVisualShowSourceSize (sources); + psphotVisualShowSourceSize (readout, sources); return true; @@ -319,2 +259,108 @@ return nSigma; } + +bool psphotMaskCosmicRay_New (psImage *mask, pmSource *source, psMaskType maskVal, psMaskType crMask) { + + // replace the source flux + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + source->mode &= ~PM_SOURCE_MODE_SUBTRACTED; + + // flag this as a CR + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + pmPeak *peak = source->peak; + psAssert (peak, "NULL peak"); + + // grab the matching footprint + pmFootprint *footprint = peak->footprint; + if (!footprint) { + // if we have not footprint, use the old code to mask by isophot + psphotMaskCosmicRay_Old (source, maskVal, crMask); + return true; + } + + if (!footprint->spans) { + // if we have not footprint, use the old code to mask by isophot + psphotMaskCosmicRay_Old (source, maskVal, crMask); + return true; + } + + // mask all of the pixels covered by the spans of the footprint + for (int j = 1; j < footprint->spans->n; j++) { + pmSpan *span1 = footprint->spans->data[j]; + + int iy = span1->y; + int xs = span1->x0; + int xe = span1->x1; + + for (int ix = xs; ix < xe; ix++) { + mask->data.U8[iy][ix] |= crMask; + } + } + return true; +} + +bool psphotMaskCosmicRay_Old (pmSource *source, psMaskType maskVal, psMaskType crMask) { + + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + pmPeak *peak = source->peak; + psAssert (peak, "NULL peak"); + + psImage *mask = source->maskView; + psImage *pixels = source->pixels; + psImage *weight = source->weight; + + // XXX This should be a recipe variable +# define SN_LIMIT 5.0 + + int xo = peak->x - pixels->col0; + int yo = peak->y - pixels->row0; + + // mark the pixels in this row to the left, then the right + for (int ix = xo; ix >= 0; ix--) { + float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]); + if (SN > SN_LIMIT) { + mask->data.U8[yo][ix] |= crMask; + } + } + for (int ix = xo + 1; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]); + if (SN > SN_LIMIT) { + mask->data.U8[yo][ix] |= crMask; + } + } + + // for each of the neighboring rows, mark the high pixels if they have a marked neighbor + // first go up: + for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) { + // mark the pixels in this row to the left, then the right + for (int ix = 0; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]); + if (SN < SN_LIMIT) continue; + + bool valid = false; + valid |= (mask->data.U8[iy+1][ix] & crMask); + valid |= (ix > 0) ? (mask->data.U8[iy+1][ix-1] & crMask) : 0; + valid |= (ix <= mask->numCols) ? (mask->data.U8[iy+1][ix+1] & crMask) : 0; + + if (!valid) continue; + mask->data.U8[iy][ix] |= crMask; + } + } + // next go down: + for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) { + // mark the pixels in this row to the left, then the right + for (int ix = 0; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]); + if (SN < SN_LIMIT) continue; + + bool valid = false; + valid |= (mask->data.U8[iy-1][ix] & crMask); + valid |= (ix > 0) ? (mask->data.U8[iy-1][ix-1] & crMask) : 0; + valid |= (ix <= mask->numCols) ? (mask->data.U8[iy-1][ix+1] & crMask) : 0; + + if (!valid) continue; + mask->data.U8[iy][ix] |= crMask; + } + } + return true; +} Index: trunk/psphot/src/psphotVisual.c =================================================================== --- trunk/psphot/src/psphotVisual.c (revision 20878) +++ trunk/psphot/src/psphotVisual.c (revision 20938) @@ -26,4 +26,50 @@ isVisual = mode; + return true; +} + +bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel) { + + KiiImage image; + KapaImageData data; + Coords coords; + + strcpy (coords.ctype, "RA---TAN"); + + psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); + if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) { + fprintf (stderr, "failed to get background values\n"); + return false; + } + + image.Nx = inImage->numCols; + image.Ny = inImage->numRows; + + ALLOCATE (image.data2d, float *, image.Ny); + for (int iy = 0; iy < image.Ny; iy++) { + ALLOCATE (image.data2d[iy], float, image.Nx); + for (int ix = 0; ix < image.Nx; ix++) { + image.data2d[iy][ix] = inImage->data.U8[iy][ix]; + } + } + + strcpy (data.name, name); + strcpy (data.file, name); + data.zero = -1; + data.range = 32; + data.logflux = 0; + + KiiSetChannel (kapaFD, channel); + KiiNewPicture2D (kapaFD, &image, &data, &coords); + + for (int iy = 0; iy < image.Ny; iy++) { + free (image.data2d[iy]); + } + free (image.data2d); + + psFree (stats); + psFree (rng); + return true; } @@ -1336,5 +1382,5 @@ } -bool psphotVisualShowSourceSize (psArray *sources) { +bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources) { int Noverlay, NOVERLAY; @@ -1398,4 +1444,6 @@ KiiLoadOverlay (kapa, overlay, Noverlay, "blue"); FREE (overlay); + + psphotVisualShowMask (kapa, readout->mask, "mask", 2); // pause and wait for user input: