Index: /branches/eam_branch_20081214/psphot/src/psphotSourceSize.c =================================================================== --- /branches/eam_branch_20081214/psphot/src/psphotSourceSize.c (revision 20977) +++ /branches/eam_branch_20081214/psphot/src/psphotSourceSize.c (revision 20978) @@ -2,8 +2,11 @@ # include +// XXX This should be a recipe variable +# define SN_LIMIT 5.0 + 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_Old (pmReadout *readout, pmSource *source, psMaskType maskVal, psMaskType crMask); bool psphotMaskCosmicRay_New (psImage *mask, pmSource *source, psMaskType maskVal, psMaskType crMask); @@ -44,4 +47,6 @@ soft = 0.0; } + + psphotVisualShowSNLimit (readout, SN_LIMIT); // loop over all source @@ -161,5 +166,6 @@ if (source->crNsigma > CR_NSIGMA_LIMIT) { // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask); - psphotMaskCosmicRay_Old (source, maskVal, crMask); + // psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask); + psphotMaskCosmicRay_Old (readout, source, maskVal, crMask); } } @@ -275,5 +281,6 @@ if (!footprint) { // if we have not footprint, use the old code to mask by isophot - psphotMaskCosmicRay_Old (source, maskVal, crMask); + // XXX this is broken, but we probably will not use it + psphotMaskCosmicRay_Old (NULL, source, maskVal, crMask); return true; } @@ -281,5 +288,6 @@ if (!footprint->spans) { // if we have not footprint, use the old code to mask by isophot - psphotMaskCosmicRay_Old (source, maskVal, crMask); + // XXX this is broken, but we probably will not use it + psphotMaskCosmicRay_Old (NULL, source, maskVal, crMask); return true; } @@ -300,5 +308,5 @@ } -bool psphotMaskCosmicRay_Old (pmSource *source, psMaskType maskVal, psMaskType crMask) { +bool psphotMaskCosmicRay_Old (pmReadout *readout, pmSource *source, psMaskType maskVal, psMaskType crMask) { source->mode |= PM_SOURCE_MODE_CR_LIMIT; @@ -306,61 +314,166 @@ 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 + // XXX should we be doing this on the smoothed image?? + psImage *pixels = readout->image; + psImage *mask = readout->mask; + psImage *weight = readout->weight; int xo = peak->x - pixels->col0; int yo = peak->y - pixels->row0; - - // mark the pixels in this row to the left, then the right + bool rowDone; + + psVector *xf = psVectorAllocEmpty (64, PS_TYPE_F32); + psVector *yf = psVectorAllocEmpty (64, PS_TYPE_F32); + + // mark the pixels in this row to the left, then the right. record the bounds of the valid pixel (SN >= SN_LIMIT) + int xs = 0; + int xe = pixels->numCols; 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; - } + if (SN < SN_LIMIT) { + xs = ix + 1; + break; + } + mask->data.U8[yo][ix] |= crMask; + psVectorAppend (xf, (float) ix); + psVectorAppend (yf, (float) yo); } 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; - } - } + if (SN < SN_LIMIT) { + xe = ix; + break; + } + mask->data.U8[yo][ix] |= crMask; + psVectorAppend (xf, (float) ix); + psVectorAppend (yf, (float) yo); + } + psphotVisualShowCRMask (xf, yf); + + int startXs = xs; + int startXe = xe; + rowDone = false; // 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; - } - } + // first go down: + for (int iy = yo - 1; !rowDone && (iy >= 0); iy--) { + rowDone = true; // we are not done if any valid pixels are found a row + + // given the range xs to xe (valid pixel range from previous row), find the valid + // pixels within this range, and extend if the ends are valid + + // newXs,newXe will define the valid pixel range for this row + int newXs = xe; + int newXe = xs; + + // first, try all of the pixels in the valid range (xs < x < xe guaranteed to be real pixels) + for (int ix = xs; ix < xe; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]); + if (SN < SN_LIMIT) { + continue; + } + newXs = PS_MIN (newXs, ix); + newXe = PS_MAX (newXe, ix + 1); + rowDone = false; + mask->data.U8[iy][ix] |= crMask; + psVectorAppend (xf, (float) ix); + psVectorAppend (yf, (float) iy); + } + + // next, try the pixels to the left of the valid range for the previous row + for (int ix = xs - 1; ix >= 0; ix--) { + float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]); + if (SN < SN_LIMIT) { + newXs = ix + 1; + break; + } + rowDone = false; + mask->data.U8[iy][ix] |= crMask; + psVectorAppend (xf, (float) ix); + psVectorAppend (yf, (float) iy); + } + + // finally, try the pixels to the right of the valid range for the previous row + for (int ix = xe; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]); + if (SN < SN_LIMIT) { + newXe = ix; + break; + } + rowDone = false; + mask->data.U8[iy][ix] |= crMask; + psVectorAppend (xf, (float) ix); + psVectorAppend (yf, (float) iy); + } + + xs = newXs; + xe = newXe; + } + psphotVisualShowCRMask (xf, yf); + + xs = startXs; + xe = startXe; + rowDone = false; + + // for each of the neighboring rows, mark the high pixels if they have a marked neighbor + // first go down: + for (int iy = yo + 1; !rowDone && (iy < pixels->numRows); iy++) { + rowDone = true; // we are not done if any valid pixels are found a row + + // given the range xs to xe (valid pixel range from previous row), find the valid + // pixels within this range, and extend if the ends are valid + + // newXs,newXe will define the valid pixel range for this row + int newXs = xe; + int newXe = xs; + + // first, try all of the pixels in the valid range (xs < x < xe guaranteed to be real pixels) + for (int ix = xs; ix < xe; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]); + if (SN < SN_LIMIT) { + continue; + } + newXs = PS_MIN (newXs, ix); + newXe = PS_MAX (newXe, ix + 1); + rowDone = false; + mask->data.U8[iy][ix] |= crMask; + psVectorAppend (xf, (float) ix); + psVectorAppend (yf, (float) iy); + } + + // next, try the pixels to the left of the valid range for the previous row + for (int ix = xs - 1; ix >= 0; ix--) { + float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]); + if (SN < SN_LIMIT) { + newXs = ix + 1; + break; + } + rowDone = false; + mask->data.U8[iy][ix] |= crMask; + psVectorAppend (xf, (float) ix); + psVectorAppend (yf, (float) iy); + } + + // finally, try the pixels to the right of the valid range for the previous row + for (int ix = xe; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]); + if (SN < SN_LIMIT) { + newXe = ix; + break; + } + rowDone = false; + mask->data.U8[iy][ix] |= crMask; + psVectorAppend (xf, (float) ix); + psVectorAppend (yf, (float) iy); + } + + xs = newXs; + xe = newXe; + } + + psphotVisualShowCRMask (xf, yf); + psFree (xf); + psFree (yf); + return true; }