Index: /branches/eam_branches/20091201/psModules/src/camera/pmFPA.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/camera/pmFPA.c (revision 26532) +++ /branches/eam_branches/20091201/psModules/src/camera/pmFPA.c (revision 26533) @@ -105,5 +105,5 @@ psFree(fpa->concepts); psFree(fpa->analysis); - psFree((void *)fpa->camera); + psFree(fpa->camera); psFree(fpa->hdu); Index: /branches/eam_branches/20091201/psModules/src/objects/pmFootprintCullPeaks.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/objects/pmFootprintCullPeaks.c (revision 26532) +++ /branches/eam_branches/20091201/psModules/src/objects/pmFootprintCullPeaks.c (revision 26533) @@ -29,11 +29,11 @@ */ -# define IN_PEAK 1 +# define IN_PEAK 1 psErrorCode pmFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint - const psImage *weight, // corresponding variance image - pmFootprint *fp, // Footprint containing mortal peaks - const float nsigma_delta, // how many sigma above local background a peak needs to be to survive - const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance - const float min_threshold) { // minimum permitted coll height + const psImage *weight, // corresponding variance image + pmFootprint *fp, // Footprint containing mortal peaks + const float nsigma_delta, // how many sigma above local background a peak needs to be to survive + const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance + const float min_threshold) { // minimum permitted coll height assert (img != NULL); assert (img->type.type == PS_TYPE_F32); assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32); @@ -42,8 +42,8 @@ if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do - return PS_ERR_NONE; - } - - psRegion subRegion; // desired subregion; 1 larger than bounding box (grr) + return PS_ERR_NONE; + } + + psRegion subRegion; // desired subregion; 1 larger than bounding box (grr) subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1; subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1; @@ -55,5 +55,5 @@ psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32); - // We need a psArray of peaks brighter than the current peak. + // We need a psArray of peaks brighter than the current peak. // We reject peaks which either: // 1) are below the local threshold @@ -67,68 +67,68 @@ // The brightest peak is always safe; go through other peaks trying to cull them for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop - const pmPeak *peak = fp->peaks->data[i]; - int x = peak->x - subImg->col0; - int y = peak->y - subImg->row0; - // - // Find the level nsigma below the peak that must separate the peak - // from any of its friends - // - assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows); - - // const float stdev = sqrt(subWt->data.F32[y][x]); - // float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev; - - const float stdev = sqrt(subWt->data.F32[y][x]); - const float flux = subImg->data.F32[y][x]; - const float fStdev = fabs(stdev/flux); - const float stdev_pad = fabs(flux * hypot(fStdev, fPadding)); - // if flux is negative, careful with fStdev - - float threshold = flux - nsigma_delta*stdev_pad; - - if (isnan(threshold) || threshold < min_threshold) { - // min_threshold is assumed to be below the detection threshold, - // so all the peaks are pmFootprint, and this isn't the brightest - continue; - } - - // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true? - if (threshold > subImg->data.F32[y][x]) { - threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON; - } - - // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint - // perhaps this should alloc a single ID image above and pass it in to be set. - - // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31) - - // XXX this is not quite there yet: - // pmFootprint *peakFootprint = NULL; - // int area = subImg->numCols * subImg->numRows; - // if (area > 30000) { - - pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x); - - // at this point brightPeaks only has the peaks brighter than the current - - // XXX need to supply the image here - // we set the IDs to either 1 (in peak) or 0 (not in peak) - pmSetFootprintID (idImg, peakFootprint, IN_PEAK); - psFree(peakFootprint); - - // Check if any of the previous (brighter) peaks are within the footprint of this peak - // If so, the current peak is bogus; drop it. - bool keep = true; - for (int j = 0; keep && (j < brightPeaks->n); j++) { - const pmPeak *peak2 = fp->peaks->data[j]; - int x2 = peak2->x - subImg->col0; - int y2 = peak2->y - subImg->row0; - if (idImg->data.S32[y2][x2] == IN_PEAK) - // There's a brighter peak within the footprint above threshold; so cull our initial peak - keep = false; - } - if (!keep) continue; - - psArrayAdd (brightPeaks, 128, fp->peaks->data[i]); + const pmPeak *peak = fp->peaks->data[i]; + int x = peak->x - subImg->col0; + int y = peak->y - subImg->row0; + // + // Find the level nsigma below the peak that must separate the peak + // from any of its friends + // + assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows); + + // const float stdev = sqrt(subWt->data.F32[y][x]); + // float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev; + + const float stdev = sqrt(subWt->data.F32[y][x]); + const float flux = subImg->data.F32[y][x]; + const float fStdev = fabs(stdev/flux); + const float stdev_pad = fabs(flux * hypot(fStdev, fPadding)); + // if flux is negative, careful with fStdev + + float threshold = flux - nsigma_delta*stdev_pad; + + if (isnan(threshold) || threshold < min_threshold) { + // min_threshold is assumed to be below the detection threshold, + // so all the peaks are pmFootprint, and this isn't the brightest + continue; + } + + // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true? + if (threshold > subImg->data.F32[y][x]) { + threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON; + } + + // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint + // perhaps this should alloc a single ID image above and pass it in to be set. + + // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31) + + // XXX this is not quite there yet: + // pmFootprint *peakFootprint = NULL; + // int area = subImg->numCols * subImg->numRows; + // if (area > 30000) { + + pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x); + + // at this point brightPeaks only has the peaks brighter than the current + + // XXX need to supply the image here + // we set the IDs to either 1 (in peak) or 0 (not in peak) + pmSetFootprintID (idImg, peakFootprint, IN_PEAK); + psFree(peakFootprint); + + // Check if any of the previous (brighter) peaks are within the footprint of this peak + // If so, the current peak is bogus; drop it. + bool keep = true; + for (int j = 0; keep && (j < brightPeaks->n); j++) { + const pmPeak *peak2 = fp->peaks->data[j]; + int x2 = peak2->x - subImg->col0; + int y2 = peak2->y - subImg->row0; + if (idImg->data.S32[y2][x2] == IN_PEAK) + // There's a brighter peak within the footprint above threshold; so cull our initial peak + keep = false; + } + if (!keep) continue; + + psArrayAdd (brightPeaks, 128, fp->peaks->data[i]); } @@ -151,9 +151,9 @@ */ psErrorCode pmFootprintCullPeaks_OLD(const psImage *img, // the image wherein lives the footprint - const psImage *weight, // corresponding variance image - pmFootprint *fp, // Footprint containing mortal peaks - const float nsigma_delta, // how many sigma above local background a peak needs to be to survive - const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance - const float min_threshold) { // minimum permitted coll height + const psImage *weight, // corresponding variance image + pmFootprint *fp, // Footprint containing mortal peaks + const float nsigma_delta, // how many sigma above local background a peak needs to be to survive + const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance + const float min_threshold) { // minimum permitted coll height assert (img != NULL); assert (img->type.type == PS_TYPE_F32); assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32); @@ -162,8 +162,8 @@ if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do - return PS_ERR_NONE; - } - - psRegion subRegion; // desired subregion; 1 larger than bounding box (grr) + return PS_ERR_NONE; + } + + psRegion subRegion; // desired subregion; 1 larger than bounding box (grr) subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1; subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1; @@ -185,79 +185,79 @@ // for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop - const pmPeak *peak = fp->peaks->data[i]; - int x = peak->x - subImg->col0; - int y = peak->y - subImg->row0; - // - // Find the level nsigma below the peak that must separate the peak - // from any of its friends - // - assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows); - - // stdev ~ sqrt(flux) : for bright regions (f ~ 1e6, sqrt(f) = 1e3 -- unrealistically small deviations) - // add additional padding in quadrature: - - const float stdev = sqrt(subWt->data.F32[y][x]); - const float flux = subImg->data.F32[y][x]; - const float fStdev = stdev/flux; - const float stdev_pad = flux * hypot(fStdev, fPadding); - float threshold = flux - nsigma_delta*stdev_pad; - - if (isnan(threshold) || threshold < min_threshold) { -#if 1 // min_threshold is assumed to be below the detection threshold, - // so all the peaks are pmFootprint, and this isn't the brightest - // XXX mark peak to be dropped - (void)psArrayRemoveIndex(fp->peaks, i); - i--; // we moved everything down one - continue; + const pmPeak *peak = fp->peaks->data[i]; + int x = peak->x - subImg->col0; + int y = peak->y - subImg->row0; + // + // Find the level nsigma below the peak that must separate the peak + // from any of its friends + // + assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows); + + // stdev ~ sqrt(flux) : for bright regions (f ~ 1e6, sqrt(f) = 1e3 -- unrealistically small deviations) + // add additional padding in quadrature: + + const float stdev = sqrt(subWt->data.F32[y][x]); + const float flux = subImg->data.F32[y][x]; + const float fStdev = stdev/flux; + const float stdev_pad = flux * hypot(fStdev, fPadding); + float threshold = flux - nsigma_delta*stdev_pad; + + if (isnan(threshold) || threshold < min_threshold) { +#if 1 // min_threshold is assumed to be below the detection threshold, + // so all the peaks are pmFootprint, and this isn't the brightest + // XXX mark peak to be dropped + (void)psArrayRemoveIndex(fp->peaks, i); + i--; // we moved everything down one + continue; #else #error n.b. We will be running LOTS of checks at this threshold, so only find the footprint once - threshold = min_threshold; + threshold = min_threshold; #endif - } - - // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true? - if (threshold > subImg->data.F32[y][x]) { - threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON; - } - - // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint - // perhaps this should alloc a single ID image above and pass it in to be set. - - const int peak_id = 1; // the ID for the peak of interest - brightPeaks->n = i; // only stop at a peak brighter than we are - - // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31) - - pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x); - brightPeaks->n = 0; // don't double free - psImage *idImg = pmSetFootprintID(NULL, peakFootprint, peak_id); - psFree(peakFootprint); - - // Check if any of the previous (brighter) peaks are within the footprint of this peak - // If so, the current peak is bogus; drop it. - int j; - for (j = 0; j < i; j++) { - const pmPeak *peak2 = fp->peaks->data[j]; - int x2 = peak2->x - subImg->col0; - int y2 = peak2->y - subImg->row0; - const int peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak - - if (peak2_id == peak_id) { // There's a brighter peak within the footprint above - ; // threshold; so cull our initial peak - (void)psArrayRemoveIndex(fp->peaks, i); - i--; // we moved everything down one - break; - } - } - if (j == i) { - j++; - } - - psFree(idImg); + } + + // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true? + if (threshold > subImg->data.F32[y][x]) { + threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON; + } + + // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint + // perhaps this should alloc a single ID image above and pass it in to be set. + + const int peak_id = 1; // the ID for the peak of interest + brightPeaks->n = i; // only stop at a peak brighter than we are + + // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31) + + pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x); + brightPeaks->n = 0; // don't double free + psImage *idImg = pmSetFootprintID(NULL, peakFootprint, peak_id); + psFree(peakFootprint); + + // Check if any of the previous (brighter) peaks are within the footprint of this peak + // If so, the current peak is bogus; drop it. + int j; + for (j = 0; j < i; j++) { + const pmPeak *peak2 = fp->peaks->data[j]; + int x2 = peak2->x - subImg->col0; + int y2 = peak2->y - subImg->row0; + const int peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak + + if (peak2_id == peak_id) { // There's a brighter peak within the footprint above + ; // threshold; so cull our initial peak + (void)psArrayRemoveIndex(fp->peaks, i); + i--; // we moved everything down one + break; + } + } + if (j == i) { + j++; + } + + psFree(idImg); } brightPeaks->n = 0; psFree(brightPeaks); - psFree((psImage *)subImg); - psFree((psImage *)subWt); + psFree(subImg); + psFree(subWt); return PS_ERR_NONE; Index: /branches/eam_branches/20091201/psModules/src/objects/pmFootprintFindAtPoint.c =================================================================== --- /branches/eam_branches/20091201/psModules/src/objects/pmFootprintFindAtPoint.c (revision 26532) +++ /branches/eam_branches/20091201/psModules/src/objects/pmFootprintFindAtPoint.c (revision 26533) @@ -29,5 +29,5 @@ * so we set appropriate mask bits * - * EAM : these function were confusingly using "startspan" and "spartspan" + * EAM : these function were confusingly using "startspan" and "spartspan" * I've rationalized them all to 'startspan' */ @@ -36,18 +36,18 @@ // An enum for what we should do with a Startspan // -typedef enum {PM_SSPAN_DOWN = 0, // scan down from this span - PM_SSPAN_UP, // scan up from this span - PM_SSPAN_RESTART, // restart scanning from this span - PM_SSPAN_DONE // this span is processed -} PM_SSPAN_DIR; // How to continue searching +typedef enum {PM_SSPAN_DOWN = 0, // scan down from this span + PM_SSPAN_UP, // scan up from this span + PM_SSPAN_RESTART, // restart scanning from this span + PM_SSPAN_DONE // this span is processed +} PM_SSPAN_DIR; // How to continue searching // // An enum for mask's pixel values. We're looking for pixels that are above threshold, and // we keep extra book-keeping information in the PM_SSPAN_STOP plane. It's simpler to be -// able to check for +// able to check for // enum { - PM_SSPAN_INITIAL = 0x0, // initial state of pixels. - PM_SSPAN_DETECTED = 0x1, // we've seen this pixel - PM_SSPAN_STOP = 0x2 // you may stop searching when you see this pixel + PM_SSPAN_INITIAL = 0x0, // initial state of pixels. + PM_SSPAN_DETECTED = 0x1, // we've seen this pixel + PM_SSPAN_STOP = 0x2 // you may stop searching when you see this pixel }; // @@ -55,17 +55,17 @@ // typedef struct { - const pmSpan *span; // save the pixel range - PM_SSPAN_DIR direction; // How to continue searching - bool stop; // should we stop searching? + const pmSpan *span; // save the pixel range + PM_SSPAN_DIR direction; // How to continue searching + bool stop; // should we stop searching? } Startspan; static void startspanFree(Startspan *sspan) { - psFree((void *)sspan->span); + psFree(sspan->span); } static Startspan * -StartspanAlloc(const pmSpan *span, // The span in question - psImage *mask, // Pixels that we've already detected - const PM_SSPAN_DIR dir // Should we continue searching towards the top of the image? +StartspanAlloc(const pmSpan *span, // The span in question + psImage *mask, // Pixels that we've already detected + const PM_SSPAN_DIR dir // Should we continue searching towards the top of the image? ) { Startspan *sspan = psAlloc(sizeof(Startspan)); @@ -75,16 +75,16 @@ sspan->direction = dir; sspan->stop = false; - - if (mask != NULL) { // remember that we've detected these pixels - psImageMaskType *mpix = &mask->data.PS_TYPE_IMAGE_MASK_DATA[span->y - mask->row0][span->x0 - mask->col0]; - - for (int i = 0; i <= span->x1 - span->x0; i++) { - mpix[i] |= PM_SSPAN_DETECTED; - if (mpix[i] & PM_SSPAN_STOP) { - sspan->stop = true; - } - } - } - + + if (mask != NULL) { // remember that we've detected these pixels + psImageMaskType *mpix = &mask->data.PS_TYPE_IMAGE_MASK_DATA[span->y - mask->row0][span->x0 - mask->col0]; + + for (int i = 0; i <= span->x1 - span->x0; i++) { + mpix[i] |= PM_SSPAN_DETECTED; + if (mpix[i] & PM_SSPAN_STOP) { + sspan->stop = true; + } + } + } + return sspan; } @@ -94,22 +94,22 @@ // static bool add_startspan(psArray *startspans, // the saved Startspans - const pmSpan *sp, // the span in question - psImage *mask, // mask of detected/stop pixels - const PM_SSPAN_DIR dir) { // the desired direction to search + const pmSpan *sp, // the span in question + psImage *mask, // mask of detected/stop pixels + const PM_SSPAN_DIR dir) { // the desired direction to search if (dir == PM_SSPAN_RESTART) { - if (add_startspan(startspans, sp, mask, PM_SSPAN_UP) || - add_startspan(startspans, sp, NULL, PM_SSPAN_DOWN)) { - return true; - } + if (add_startspan(startspans, sp, mask, PM_SSPAN_UP) || + add_startspan(startspans, sp, NULL, PM_SSPAN_DOWN)) { + return true; + } } else { - Startspan *sspan = StartspanAlloc(sp, mask, dir); - if (sspan->stop) { // we detected a stop bit - psFree(sspan); // don't allocate new span - - return true; - } else { - psArrayAdd(startspans, 1, sspan); - psFree(sspan); // as it's now owned by startspans - } + Startspan *sspan = StartspanAlloc(sp, mask, dir); + if (sspan->stop) { // we detected a stop bit + psFree(sspan); // don't allocate new span + + return true; + } else { + psArrayAdd(startspans, 1, sspan); + psFree(sspan); // as it's now owned by startspans + } } @@ -127,44 +127,44 @@ */ static bool do_startspan(pmFootprint *fp, // the footprint that we're building - const psImage *img, // the psImage we're working on - psImage *mask, // the associated masks - const float threshold, // Threshold - psArray *startspans) { // specify which span to process next - bool F32 = false; // is this an F32 image? + const psImage *img, // the psImage we're working on + psImage *mask, // the associated masks + const float threshold, // Threshold + psArray *startspans) { // specify which span to process next + bool F32 = false; // is this an F32 image? if (img->type.type == PS_TYPE_F32) { - F32 = true; + F32 = true; } else if (img->type.type == PS_TYPE_S32) { - F32 = false; - } else { // N.b. You can't trivially add more cases here; F32 is just a bool - psError(PS_ERR_UNKNOWN, true, "Unsupported psImage type: %d", img->type.type); - return NULL; - } - - psF32 *imgRowF32 = NULL; // row pointer if F32 - psS32 *imgRowS32 = NULL; // " " " " !F32 - psImageMaskType *maskRow = NULL; // masks's row pointer - + F32 = false; + } else { // N.b. You can't trivially add more cases here; F32 is just a bool + psError(PS_ERR_UNKNOWN, true, "Unsupported psImage type: %d", img->type.type); + return NULL; + } + + psF32 *imgRowF32 = NULL; // row pointer if F32 + psS32 *imgRowS32 = NULL; // " " " " !F32 + psImageMaskType *maskRow = NULL; // masks's row pointer + const int row0 = img->row0; const int col0 = img->col0; const int numRows = img->numRows; const int numCols = img->numCols; - + /********************************************************************************************************/ - + Startspan *sspan = NULL; for (int i = 0; i < startspans->n; i++) { - sspan = startspans->data[i]; - if (sspan->direction != PM_SSPAN_DONE) { - break; - } - if (sspan->stop) { - break; - } + sspan = startspans->data[i]; + if (sspan->direction != PM_SSPAN_DONE) { + break; + } + if (sspan->stop) { + break; + } } if (sspan == NULL || sspan->direction == PM_SSPAN_DONE) { // no more Startspans to process - return false; - } - if (sspan->stop) { // they don't want any more spans processed - return false; + return false; + } + if (sspan->stop) { // they don't want any more spans processed + return false; } /* @@ -179,111 +179,111 @@ * Go through image identifying objects */ - int nx0, nx1 = -1; // new values of x0, x1 + int nx0, nx1 = -1; // new values of x0, x1 const int di = (dir == PM_SSPAN_UP) ? 1 : -1; // how much i changes to get to the next row - bool stop = false; // should I stop searching for spans? + bool stop = false; // should I stop searching for spans? for (int i = sspan->span->y -row0 + di; i < numRows && i >= 0; i += di) { - imgRowF32 = img->data.F32[i]; // only one of - imgRowS32 = img->data.S32[i]; // these is valid! - maskRow = mask->data.PS_TYPE_IMAGE_MASK_DATA[i]; - // - // Search left from the pixel diagonally to the left of (i - di, x0). If there's - // a connected span there it may need to grow up and/or down, so push it onto - // the stack for later consideration - // - nx0 = -1; - for (int j = x0 - 1; j >= -1; j--) { - double pixVal = (j < 0) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]); - if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) { - if (j < x0 - 1) { // we found some pixels above threshold - nx0 = j + 1; - } - break; - } - } - - if (nx0 < 0) { // no span to the left - nx1 = x0 - 1; // we're going to resume searching at nx1 + 1 - } else { - // - // Search right in leftmost span - // - //nx1 = 0; // make gcc happy - for (int j = nx0 + 1; j <= numCols; j++) { - double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]); - if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) { - nx1 = j - 1; - break; - } - } - - const pmSpan *sp = pmFootprintAddSpan(fp, i + row0, nx0 + col0, nx1 + col0); - - if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) { - stop = true; - break; - } - } - // - // Now look for spans connected to the old span. The first of these we'll - // simply process, but others will have to be deferred for later consideration. - // - // In fact, if the span overhangs to the right we'll have to defer the overhang - // until later too, as it too can grow in both directions - // - // Note that column numCols exists virtually, and always ends the last span; this - // is why we claim below that sx1 is always set - // - bool first = false; // is this the first new span detected? - for (int j = nx1 + 1; j <= x1 + 1; j++) { - double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]); - if (!(maskRow[j] & PM_SSPAN_DETECTED) && pixVal >= threshold) { - int sx0 = j++; // span that we're working on is sx0:sx1 - int sx1 = -1; // We know that if we got here, we'll also set sx1 - for (; j <= numCols; j++) { - double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]); - if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) { // end of span - sx1 = j; - break; - } - } - assert (sx1 >= 0); - - const pmSpan *sp; - if (first) { - if (sx1 <= x1) { - sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1); - if (add_startspan(startspans, sp, mask, PM_SSPAN_DONE)) { - stop = true; - break; - } - } else { // overhangs to right - sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, x1 + col0); - if (add_startspan(startspans, sp, mask, PM_SSPAN_DONE)) { - stop = true; - break; - } - sp = pmFootprintAddSpan(fp, i + row0, x1 + 1 + col0, sx1 + col0 - 1); - if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) { - stop = true; - break; - } - } - first = false; - } else { - sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1); - if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) { - stop = true; - break; - } - } - } - } - - if (stop || first == false) { // we're done - break; - } - - x0 = nx0; x1 = nx1; + imgRowF32 = img->data.F32[i]; // only one of + imgRowS32 = img->data.S32[i]; // these is valid! + maskRow = mask->data.PS_TYPE_IMAGE_MASK_DATA[i]; + // + // Search left from the pixel diagonally to the left of (i - di, x0). If there's + // a connected span there it may need to grow up and/or down, so push it onto + // the stack for later consideration + // + nx0 = -1; + for (int j = x0 - 1; j >= -1; j--) { + double pixVal = (j < 0) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]); + if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) { + if (j < x0 - 1) { // we found some pixels above threshold + nx0 = j + 1; + } + break; + } + } + + if (nx0 < 0) { // no span to the left + nx1 = x0 - 1; // we're going to resume searching at nx1 + 1 + } else { + // + // Search right in leftmost span + // + //nx1 = 0; // make gcc happy + for (int j = nx0 + 1; j <= numCols; j++) { + double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]); + if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) { + nx1 = j - 1; + break; + } + } + + const pmSpan *sp = pmFootprintAddSpan(fp, i + row0, nx0 + col0, nx1 + col0); + + if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) { + stop = true; + break; + } + } + // + // Now look for spans connected to the old span. The first of these we'll + // simply process, but others will have to be deferred for later consideration. + // + // In fact, if the span overhangs to the right we'll have to defer the overhang + // until later too, as it too can grow in both directions + // + // Note that column numCols exists virtually, and always ends the last span; this + // is why we claim below that sx1 is always set + // + bool first = false; // is this the first new span detected? + for (int j = nx1 + 1; j <= x1 + 1; j++) { + double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]); + if (!(maskRow[j] & PM_SSPAN_DETECTED) && pixVal >= threshold) { + int sx0 = j++; // span that we're working on is sx0:sx1 + int sx1 = -1; // We know that if we got here, we'll also set sx1 + for (; j <= numCols; j++) { + double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]); + if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) { // end of span + sx1 = j; + break; + } + } + assert (sx1 >= 0); + + const pmSpan *sp; + if (first) { + if (sx1 <= x1) { + sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1); + if (add_startspan(startspans, sp, mask, PM_SSPAN_DONE)) { + stop = true; + break; + } + } else { // overhangs to right + sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, x1 + col0); + if (add_startspan(startspans, sp, mask, PM_SSPAN_DONE)) { + stop = true; + break; + } + sp = pmFootprintAddSpan(fp, i + row0, x1 + 1 + col0, sx1 + col0 - 1); + if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) { + stop = true; + break; + } + } + first = false; + } else { + sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1); + if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) { + stop = true; + break; + } + } + } + } + + if (stop || first == false) { // we're done + break; + } + + x0 = nx0; x1 = nx1; } /* @@ -309,22 +309,22 @@ */ pmFootprint * -pmFootprintsFindAtPoint(const psImage *img, // image to search - const float threshold, // Threshold - const psArray *peaks, // array of peaks; finding one terminates search for footprint - int row, int col) { // starting position (in img's parent's coordinate system) +pmFootprintsFindAtPoint(const psImage *img, // image to search + const float threshold, // Threshold + const psArray *peaks, // array of peaks; finding one terminates search for footprint + int row, int col) { // starting position (in img's parent's coordinate system) assert(img != NULL); - bool F32 = false; // is this an F32 image? + bool F32 = false; // is this an F32 image? if (img->type.type == PS_TYPE_F32) { F32 = true; } else if (img->type.type == PS_TYPE_S32) { F32 = false; - } else { // N.b. You can't trivially add more cases here; F32 is just a bool + } else { // N.b. You can't trivially add more cases here; F32 is just a bool psError(PS_ERR_UNKNOWN, true, "Unsupported psImage type: %d", img->type.type); return NULL; } - psF32 *imgRowF32 = NULL; // row pointer if F32 - psS32 *imgRowS32 = NULL; // " " " " !F32 - + psF32 *imgRowF32 = NULL; // row pointer if F32 + psS32 *imgRowS32 = NULL; // " " " " !F32 + const int row0 = img->row0; const int col0 = img->col0; @@ -339,5 +339,5 @@ psError(PS_ERR_BAD_PARAMETER_VALUE, true, "row/col == (%d, %d) are out of bounds [%d--%d, %d--%d]", - row + row0, col + col0, row0, row0 + numRows - 1, col0, col0 + numCols - 1); + row + row0, col + col0, row0, row0 + numRows - 1, col0, col0 + numCols - 1); return NULL; } @@ -347,5 +347,5 @@ return pmFootprintAlloc(0, img); } - + pmFootprint *fp = pmFootprintAlloc(1 + img->numRows/10, img); /* @@ -364,6 +364,6 @@ if (peaks != NULL) { for (int i = 0; i < peaks->n; i++) { - pmPeak *peak = peaks->data[i]; - mask->data.PS_TYPE_IMAGE_MASK_DATA[peak->y - mask->row0][peak->x - mask->col0] |= PM_SSPAN_STOP; + pmPeak *peak = peaks->data[i]; + mask->data.PS_TYPE_IMAGE_MASK_DATA[peak->y - mask->row0][peak->x - mask->col0] |= PM_SSPAN_STOP; } } @@ -372,22 +372,22 @@ */ psArray *startspans = psArrayAllocEmpty(1); // spans where we have to restart the search - - imgRowF32 = img->data.F32[row]; // only one of - imgRowS32 = img->data.S32[row]; // these is valid! + + imgRowF32 = img->data.F32[row]; // only one of + imgRowS32 = img->data.S32[row]; // these is valid! psImageMaskType *maskRow = mask->data.PS_TYPE_IMAGE_MASK_DATA[row]; { int i; for (i = col; i >= 0; i--) { - pixVal = F32 ? imgRowF32[i] : imgRowS32[i]; - if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) { - break; - } + pixVal = F32 ? imgRowF32[i] : imgRowS32[i]; + if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) { + break; + } } int i0 = i; for (i = col; i < numCols; i++) { - pixVal = F32 ? imgRowF32[i] : imgRowS32[i]; - if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) { - break; - } + pixVal = F32 ? imgRowF32[i] : imgRowS32[i]; + if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) { + break; + } } int i1 = i; @@ -404,6 +404,6 @@ */ psFree(mask); - psFree(startspans); // restores the image pixel - - return fp; // pmFootprint really + psFree(startspans); // restores the image pixel + + return fp; // pmFootprint really }