Index: trunk/psModules/src/objects/Makefile.am =================================================================== --- trunk/psModules/src/objects/Makefile.am (revision 27657) +++ trunk/psModules/src/objects/Makefile.am (revision 27672) @@ -12,4 +12,5 @@ pmFootprintCullPeaks.c \ pmFootprintFind.c \ + pmFootprintSpans.c \ pmFootprintFindAtPoint.c \ pmFootprintIDs.c \ @@ -75,4 +76,5 @@ pmSpan.h \ pmFootprint.h \ + pmFootprintSpans.h \ pmPeaks.h \ pmMoments.h \ Index: trunk/psModules/src/objects/pmFootprint.c =================================================================== --- trunk/psModules/src/objects/pmFootprint.c (revision 27657) +++ trunk/psModules/src/objects/pmFootprint.c (revision 27672) @@ -53,4 +53,5 @@ assert(nspan >= 0); footprint->npix = 0; + footprint->nspans = 0; // we may allocate more spans than we set -- this is the number of active spans footprint->spans = psArrayAllocEmpty(nspan); footprint->peaks = psArrayAlloc(0); @@ -73,4 +74,23 @@ } +bool pmFootprintAllocEmptySpans (pmFootprint *footprint, int nSpans) { + + psArrayRealloc (footprint->spans, nSpans); + for (int i = 0; i < nSpans; i++) { + footprint->spans->data[i] = pmSpanAlloc(-1, -1, -1); + } + footprint->spans->n = nSpans; + return true; +} + +// reset the footprint containers +bool pmFootprintInit(pmFootprint *footprint) { + + footprint->bbox.x0 = footprint->bbox.y0 = 0; + footprint->bbox.x1 = footprint->bbox.y1 = -1; + footprint->nspans = 0; + return true; +} + bool pmFootprintTest(const psPtr ptr) { return (psMemGetDeallocator(ptr) == (psFreeFunc)footprintFree); @@ -103,8 +123,10 @@ psArrayAdd(fp->spans, 1, sp); psFree(sp); - + + fp->nspans ++; + fp->npix += x1 - x0 + 1; - if (fp->spans->n == 1) { + if (fp->nspans == 1) { fp->bbox.x0 = x0; fp->bbox.x1 = x1; @@ -121,7 +143,55 @@ } + +// Set the next available elements of the nSpan entry in footprint->spans +pmSpan *pmFootprintSetSpan(pmFootprint *fp, // the footprint to add to + const int y, // row to add + int x0, // range of + int x1) { // columns + + if (x1 < x0) { + int tmp = x0; + x0 = x1; + x1 = tmp; + } + + int N = fp->nspans; + if (N == fp->spans->n) { + // if we need more space, extend fp->spans as needed + int Nalloc = fp->spans->n + 100; + psArrayRealloc(fp->spans, Nalloc); + fp->spans->n = Nalloc; + for (int i = N; i < fp->spans->n; i++) { + fp->spans->data[i] = pmSpanAlloc(-1, -1, -1); + } + } + + pmSpan *span = fp->spans->data[N]; + span->y = y; + span->x0 = x0; + span->x1 = x1; + + fp->nspans ++; + + fp->npix += x1 - x0 + 1; + + if (fp->nspans == 1) { + fp->bbox.x0 = x0; + fp->bbox.x1 = x1; + fp->bbox.y0 = y; + fp->bbox.y1 = y; + } else { + if (x0 < fp->bbox.x0) fp->bbox.x0 = x0; + if (x1 > fp->bbox.x1) fp->bbox.x1 = x1; + if (y < fp->bbox.y0) fp->bbox.y0 = y; + if (y > fp->bbox.y1) fp->bbox.y1 = y; + } + + return span; +} + void pmFootprintSetBBox(pmFootprint *fp) { assert (fp != NULL); - if (fp->spans->n == 0) { + if (fp->nspans == 0) { return; } @@ -132,5 +202,5 @@ int y1 = sp->y; - for (int i = 1; i < fp->spans->n; i++) { + for (int i = 1; i < fp->nspans; i++) { sp = fp->spans->data[i]; @@ -150,5 +220,5 @@ assert (fp != NULL); int npix = 0; - for (int i = 0; i < fp->spans->n; i++) { + for (int i = 0; i < fp->nspans; i++) { pmSpan *span = fp->spans->data[i]; npix += span->x1 - span->x0 + 1; Index: trunk/psModules/src/objects/pmFootprint.h =================================================================== --- trunk/psModules/src/objects/pmFootprint.h (revision 27657) +++ trunk/psModules/src/objects/pmFootprint.h (revision 27672) @@ -13,10 +13,11 @@ #include #include "pmSpan.h" - +#include "pmFootprintSpans.h" typedef struct { const int id; //!< unique ID int npix; //!< number of pixels in this pmFootprint - psArray *spans; //!< the pmSpans + int nspans; + psArray *spans; //!< the allocated pmSpans psRegion bbox; //!< the pmFootprint's bounding box psArray *peaks; //!< the peaks lying in this footprint @@ -26,5 +27,8 @@ pmFootprint *pmFootprintAlloc(int nspan, const psImage *img); +bool pmFootprintInit(pmFootprint *footprint); bool pmFootprintTest(const psPtr ptr); + +bool pmFootprintAllocEmptySpans (pmFootprint *footprint, int nSpans); pmFootprint *pmFootprintNormalize(pmFootprint *fp); @@ -37,9 +41,30 @@ int x1); // columns +pmSpan *pmFootprintSetSpan(pmFootprint *fp, // the footprint to add to + const int y, // row to add + int x0, // range of + int x1); // columns + psArray *pmFootprintsFind(const psImage *img, const float threshold, const int npixMin); -pmFootprint *pmFootprintsFindAtPoint(const psImage *img, - const float threshold, - const psArray *peaks, - int row, int col); + +bool pmFootprintsFindAtPoint(pmFootprint *fp, + pmFootprintSpans *fpSpans, + const psImage *img, // image to search + psImage *mask, + const float threshold, // Threshold + const psArray *peaks, // array of peaks; finding one terminates search for footprint + int row, int col); + +// pmFootprint *pmFootprintsFindAtPoint(const psImage *img, +// const float threshold, +// const psArray *peaks, +// int row, int col); + +bool pmFootprintSpansBuild(pmFootprint *fp, // the footprint that we're building + pmFootprintSpans *fpSpans, + const psImage *img, // the psImage we're working on + psImage *mask, // the associated masks + const float threshold // Threshold + ); psArray *pmFootprintArrayGrow(const psArray *footprints, int r); Index: trunk/psModules/src/objects/pmFootprintCullPeaks.c =================================================================== --- trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 27657) +++ trunk/psModules/src/objects/pmFootprintCullPeaks.c (revision 27672) @@ -20,5 +20,8 @@ #include "pmSpan.h" #include "pmFootprint.h" +#include "pmFootprintSpans.h" #include "pmPeaks.h" + +bool dumpfootprints (pmFootprint *fp, pmFootprintSpans *fpSp); /* @@ -65,4 +68,13 @@ psArrayAdd (brightPeaks, 128, fp->peaks->data[0]); + // allocate the peakFootprint and peakFPSpans containers -- these are re-used by pmFootprintsFindAtPoint to minimize allocs in this function + pmFootprint *peakFootprint = pmFootprintAlloc(fp->nspans, subImg); + pmFootprintSpans *peakFPSpans = pmFootprintSpansAlloc(2*fp->nspans); + + // allocate empty spans for the footprints + pmFootprintAllocEmptySpans(peakFootprint, fp->nspans); + + psImage *subMask = psImageCopy(NULL, subImg, PS_TYPE_IMAGE_MASK); + // 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 @@ -98,22 +110,16 @@ } - // 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); + // init peakFootprint here? + pmFootprintsFindAtPoint(peakFootprint, peakFPSpans, subImg, subMask, threshold, brightPeaks, peak->y, peak->x); + if (peakFPSpans->nStartSpans > 2000) { + // dumpfootprints(peakFootprint, peakFPSpans); + // fprintf (stderr, "big footprint %d : %d\n", peakFootprint->nspans, peakFPSpans->nStartSpans); + fprintf (stderr, "big test footprint: %f %f to %f %f (%d pix)\n", peakFootprint->bbox.x0, peakFootprint->bbox.y0, peakFootprint->bbox.x1, peakFootprint->bbox.y1, peakFootprint->npix); + } // 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); // If this peak has not already been assigned to a source, then we can look for any @@ -144,4 +150,7 @@ psFree(subImg); psFree(subWt); + psFree(subMask); + psFree(peakFootprint); + psFree(peakFPSpans); return PS_ERR_NONE; @@ -149,122 +158,24 @@ - /* - * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently - * isolated. More precisely, for each peak find the highest coll that you'd have to traverse - * to reach a still higher peak --- and if that coll's more than nsigma DN below your - * starting point, discard the peak. - */ -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 - assert (img != NULL); assert (img->type.type == PS_TYPE_F32); - assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32); - assert (img->row0 == weight->row0 && img->col0 == weight->col0); - assert (fp != NULL); +bool dumpfootprints (pmFootprint *fp, pmFootprintSpans *fpSp) { - if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do - return PS_ERR_NONE; + FILE *f1 = fopen ("fp.dat", "w"); + if (!f1) return false; + + for (int i = 0; i < fp->nspans; i++) { + pmSpan *sp = fp->spans->data[i]; + fprintf (f1, "%d - %d : %d\n", sp->x0, sp->x1, sp->y); } + fclose (f1); - 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; - const psImage *subImg = psImageSubset((psImage *)img, subRegion); - const psImage *subWt = psImageSubset((psImage *)weight, subRegion); - assert (subImg != NULL && subWt != NULL); - // - // We need a psArray of peaks brighter than the current peak. We'll fake this - // by reusing the fp->peaks but lying about n. - // - // We do this for efficiency (otherwise I'd need two peaks lists), and we are - // rather too chummy with psArray in consequence. But it works. - // - psArray *brightPeaks = psArrayAlloc(0); - psFree(brightPeaks->data); - brightPeaks->data = psMemIncrRefCounter(fp->peaks->data);// use the data from fp->peaks - // - // 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); + FILE *f2 = fopen ("fpSp.dat", "w"); + if (!f2) return false; - // 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; -#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); + for (int i = 0; i < fpSp->nStartSpans; i++) { + pmStartSpan *sp = fpSp->startspans->data[i]; + if (!sp->span) continue; + fprintf (f2, "%d - %d : %d\n", sp->span->x0, sp->span->x1, sp->span->y); } - - brightPeaks->n = 0; psFree(brightPeaks); - psFree(subImg); - psFree(subWt); - - return PS_ERR_NONE; + fclose (f2); + return true; } - Index: trunk/psModules/src/objects/pmFootprintFindAtPoint.c =================================================================== --- trunk/psModules/src/objects/pmFootprintFindAtPoint.c (revision 27657) +++ trunk/psModules/src/objects/pmFootprintFindAtPoint.c (revision 27672) @@ -21,103 +21,9 @@ #include "pmFootprint.h" #include "pmPeaks.h" - -/* - * A data structure to hold the starting point for a search for pixels above threshold, - * used by pmFootprintsFindAtPoint - * - * We don't want to find this span again --- it's already part of the footprint --- - * so we set appropriate mask bits - * - * EAM : these function were confusingly using "startspan" and "spartspan" - * I've rationalized them all to 'startspan' - */ - -// -// 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 -// -// 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 -// -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 -}; -// -// The struct that remembers how to [re-]start scanning the image for pixels -// -typedef struct { - 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(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? - ) { - Startspan *sspan = psAlloc(sizeof(Startspan)); - psMemSetDeallocator(sspan, (psFreeFunc)startspanFree); - - sspan->span = psMemIncrRefCounter((void *)span); - 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; - } - } - } - - return sspan; -} - -// -// Add a new Startspan to an array of Startspans. Iff we see a stop bit, return true -// -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 - if (dir == PM_SSPAN_RESTART) { - 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 - } - } - - return false; -} +#include "pmFootprintSpans.h" /************************************************************************************************************/ /* - * Search the image for pixels above threshold, starting at a single Startspan. + * Search the image for pixels above threshold, starting at a single pmStartSpan. * We search the array looking for one to process; it'd be better to move the * ones that we're done with to the end, but it probably isn't worth it for @@ -126,9 +32,10 @@ * This is the guts of pmFootprintsFindAtPoint */ -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 pmFootprintSpansBuild(pmFootprint *fp, // the footprint that we're building + pmFootprintSpans *fpSpans, + const psImage *img, // the psImage we're working on + psImage *mask, // the associated masks + const float threshold // Threshold + ) { bool F32 = false; // is this an F32 image? if (img->type.type == PS_TYPE_F32) { @@ -152,36 +59,37 @@ /********************************************************************************************************/ - Startspan *sspan = NULL; - for (int i = 0; i < startspans->n; i++) { - sspan = startspans->data[i]; - if (sspan->direction != PM_SSPAN_DONE) { + pmStartSpan *startspan = NULL; + for (int i = 0; i < fpSpans->nStartSpans; i++) { + startspan = fpSpans->startspans->data[i]; + if (startspan->direction != PM_STARTSPAN_DONE) { break; } - if (sspan->stop) { + if (startspan->stop) { break; } } - if (sspan == NULL || sspan->direction == PM_SSPAN_DONE) { // no more Startspans to process + if (startspan == NULL || startspan->direction == PM_STARTSPAN_DONE) { // no more pmStartSpans to process return false; } - if (sspan->stop) { // they don't want any more spans processed + if (startspan->stop) { // they don't want any more spans processed return false; } + /* * Work */ - const PM_SSPAN_DIR dir = sspan->direction; + const PM_STARTSPAN_DIR dir = startspan->direction; /* * Set initial span to the startspan */ - int x0 = sspan->span->x0 - col0, x1 = sspan->span->x1 - col0; + int x0 = startspan->span->x0 - col0, x1 = startspan->span->x1 - col0; /* * Go through image identifying objects */ 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 + const int di = (dir == PM_STARTSPAN_UP) ? 1 : -1; // how much i changes to get to the next row bool stop = false; // should I stop searching for spans? - for (int i = sspan->span->y -row0 + di; i < numRows && i >= 0; i += di) { + for (int i = startspan->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! @@ -195,5 +103,5 @@ 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 ((maskRow[j] & PM_STARTSPAN_DETECTED) || pixVal < threshold) { if (j < x0 - 1) { // we found some pixels above threshold nx0 = j + 1; @@ -209,8 +117,7 @@ // 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) { + if ((maskRow[j] & PM_STARTSPAN_DETECTED) || pixVal < threshold) { nx1 = j - 1; break; @@ -218,7 +125,8 @@ } - const pmSpan *sp = pmFootprintAddSpan(fp, i + row0, nx0 + col0, nx1 + col0); - - if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) { + pmSpan *sp = pmFootprintSetSpan(fp, i + row0, nx0 + col0, nx1 + col0); + bool status = pmFootprintSpansSet(fpSpans, sp, mask, PM_STARTSPAN_RESTART); + // fprintf (stderr, "set 1: %d vs %d\n", fp->nspans, fpSpans->nStartSpans); + if (status) { stop = true; break; @@ -238,10 +146,10 @@ 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) { + if (!(maskRow[j] & PM_STARTSPAN_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 + if ((maskRow[j] & PM_STARTSPAN_DETECTED) || pixVal < threshold) { // end of span sx1 = j; break; @@ -250,20 +158,26 @@ assert (sx1 >= 0); - const pmSpan *sp; + 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)) { + sp = pmFootprintSetSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1); + bool status = pmFootprintSpansSet(fpSpans, sp, mask, PM_STARTSPAN_DONE); + // fprintf (stderr, "set 2: %d vs %d\n", fp->nspans, fpSpans->nStartSpans); + if (status) { 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)) { + sp = pmFootprintSetSpan(fp, i + row0, sx0 + col0, x1 + col0); + bool status = pmFootprintSpansSet(fpSpans, sp, mask, PM_STARTSPAN_DONE); + // fprintf (stderr, "set 3: %d vs %d\n", fp->nspans, fpSpans->nStartSpans); + if (status) { stop = true; break; } - sp = pmFootprintAddSpan(fp, i + row0, x1 + 1 + col0, sx1 + col0 - 1); - if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) { + sp = pmFootprintSetSpan(fp, i + row0, x1 + 1 + col0, sx1 + col0 - 1); + status = pmFootprintSpansSet(fpSpans, sp, mask, PM_STARTSPAN_RESTART); + // fprintf (stderr, "set 4: %d vs %d\n", fp->nspans, fpSpans->nStartSpans); + if (status) { stop = true; break; @@ -272,6 +186,8 @@ first = false; } else { - sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1); - if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) { + sp = pmFootprintSetSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1); + bool status = pmFootprintSpansSet(fpSpans, sp, mask, PM_STARTSPAN_RESTART); + // fprintf (stderr, "set 5: %d vs %d\n", fp->nspans, fpSpans->nStartSpans); + if (status) { stop = true; break; @@ -291,5 +207,5 @@ */ - sspan->direction = PM_SSPAN_DONE; + startspan->direction = PM_STARTSPAN_DONE; return stop ? false : true; } @@ -307,103 +223,113 @@ * are not sorted by increasing y, x0, x1. If this matters to you, call * pmFootprintNormalize() + * + * The calling function must supply a footprint allocated with a reasonable amount of space + * */ -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) - assert(img != NULL); - - 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 - 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 - - const int row0 = img->row0; - const int col0 = img->col0; - const int numRows = img->numRows; - const int numCols = img->numCols; -/* - * Is point in image, and above threshold? - */ - row -= row0; col -= col0; - if (row < 0 || row >= numRows || - col < 0 || col >= numCols) { + +bool pmFootprintsFindAtPoint(pmFootprint *fp, + pmFootprintSpans *fpSpans, + const psImage *img, // image to search + psImage *mask, + 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) + psAssert(img, "image must be supplied"); + psAssert(fp, "footprint must be supplied"); + psAssert(fpSpans, "footprint spans must be supplied"); + + 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 + psError(PS_ERR_UNKNOWN, true, "Unsupported psImage type: %d", img->type.type); + return false; + } + psF32 *imgRowF32 = NULL; // row pointer if F32 + psS32 *imgRowS32 = NULL; // " " " " !F32 + + const int row0 = img->row0; + const int col0 = img->col0; + const int numRows = img->numRows; + const int numCols = img->numCols; + + /* + * Is point in image, and above threshold? + */ + row -= row0; col -= col0; + if (row < 0 || row >= numRows || + col < 0 || col >= numCols) { 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); - return NULL; - } - - double pixVal = F32 ? img->data.F32[row][col] : img->data.S32[row][col]; - if (pixVal < threshold) { - return pmFootprintAlloc(0, img); - } - - pmFootprint *fp = pmFootprintAlloc(1 + img->numRows/10, img); -/* - * We need a mask for two purposes; to indicate which pixels are already detected, - * and to store the "stop" pixels --- those that, once reached, should stop us - * looking for the rest of the pmFootprint. These are generally set from peaks. - */ - psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); - P_PSIMAGE_SET_ROW0(mask, row0); - P_PSIMAGE_SET_COL0(mask, col0); - psImageInit(mask, PM_SSPAN_INITIAL); - // - // Set stop bits from peaks list - // - assert (peaks == NULL || peaks->n == 0 || psMemCheckPeak(peaks->data[0])); - 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; - } - } -/* - * Find starting span passing through (row, col) - */ - 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! - 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; - } - } - int i0 = i; - for (i = col; i < numCols; i++) { - pixVal = F32 ? imgRowF32[i] : imgRowS32[i]; - if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) { - break; - } - } - int i1 = i; - const pmSpan *sp = pmFootprintAddSpan(fp, row + row0, i0 + col0 + 1, i1 + col0 - 1); - - (void)add_startspan(startspans, sp, mask, PM_SSPAN_RESTART); - } - /* - * Now workout from those Startspans, searching for pixels above threshold - */ - while (do_startspan(fp, img, mask, threshold, startspans)) continue; - /* - * Cleanup - */ - psFree(mask); - psFree(startspans); // restores the image pixel - - return fp; // pmFootprint really + return false; + } + + double pixVal = F32 ? img->data.F32[row][col] : img->data.S32[row][col]; + if (pixVal < threshold) { + return true; + } + + /* + * We need a mask for two purposes; to indicate which pixels are already detected, + * and to store the "stop" pixels --- those that, once reached, should stop us + * looking for the rest of the pmFootprint. These are generally set from peaks. + */ + + pmFootprintInit(fp); + pmFootprintSpansInit(fpSpans); + psImageInit(mask, PM_STARTSPAN_INITIAL); + + // fprintf (stderr, "init: %d vs %d\n", fp->nspans, fpSpans->nStartSpans); + + // + // Set stop bits from peaks list + // + assert (peaks == NULL || peaks->n == 0 || psMemCheckPeak(peaks->data[0])); + 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_STARTSPAN_STOP; + } + } + + /* + * Find starting span passing through (row, col) + */ + 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_STARTSPAN_DETECTED) || pixVal < threshold) { + break; + } + } + int i0 = i; + for (i = col; i < numCols; i++) { + pixVal = F32 ? imgRowF32[i] : imgRowS32[i]; + if ((maskRow[i] & PM_STARTSPAN_DETECTED) || pixVal < threshold) { + break; + } + } + int i1 = i; + pmSpan *sp = pmFootprintSetSpan(fp, row + row0, i0 + col0 + 1, i1 + col0 - 1); + pmFootprintSpansSet(fpSpans, sp, mask, PM_STARTSPAN_RESTART); + // fprintf (stderr, "first: %d vs %d\n", fp->nspans, fpSpans->nStartSpans); + } + /* + * Now workout from those pmStartSpans, searching for pixels above threshold + */ + while (pmFootprintSpansBuild(fp, fpSpans, img, mask, threshold)) continue; + /* + * Cleanup + */ + // psFree(mask); + // psFree(startspans); // restores the image pixel + + return fp; // pmFootprint really } Index: trunk/psModules/src/objects/pmFootprintIDs.c =================================================================== --- trunk/psModules/src/objects/pmFootprintIDs.c (revision 27657) +++ trunk/psModules/src/objects/pmFootprintIDs.c (revision 27672) @@ -32,5 +32,5 @@ const int row0 = fp->region.y0; - for (int j = 0; j < fp->spans->n; j++) { + for (int j = 0; j < fp->nspans; j++) { const pmSpan *span = fp->spans->data[j]; psS32 *imgRow = idImage->data.S32[span->y - row0]; Index: trunk/psModules/src/objects/pmFootprintSpans.c =================================================================== --- trunk/psModules/src/objects/pmFootprintSpans.c (revision 27672) +++ trunk/psModules/src/objects/pmFootprintSpans.c (revision 27672) @@ -0,0 +1,136 @@ +/* @file pmFootprintFindAtPoint.c + * find footprints in a small image relative to a reference point + * + * @author RHL, Princeton & IfA; EAM, IfA + * + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-01-27 06:39:38 $ + * Copyright 2006 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include "pmSpan.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmFootprintSpans.h" + +static void pmStartSpanFree(pmStartSpan *sspan) { + return; +} + +// Allocate an un-assigned pmStartSpan +pmStartSpan *pmStartSpanAlloc() { + pmStartSpan *sspan = psAlloc(sizeof(pmStartSpan)); + psMemSetDeallocator(sspan, (psFreeFunc)pmStartSpanFree); + + sspan->span = NULL; + sspan->direction = PM_STARTSPAN_NONE; + sspan->stop = false; + + return sspan; +} + +// Assign a pmSpan to this pmStartSpan +bool pmStartSpanSet(pmStartSpan *sspan, + pmSpan *span, // The span in question + psImage *mask, // Pixels that we've already detected + const PM_STARTSPAN_DIR dir // Should we continue searching towards the top of the image? + ) { + sspan->span = span; // view on the span (we do not own this memory) + sspan->direction = dir; + sspan->stop = false; + + if (mask) { // 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_STARTSPAN_DETECTED; + if (mpix[i] & PM_STARTSPAN_STOP) { + sspan->stop = true; + } + } + } + + return true; +} + +// remove a span assignment +bool pmStartSpanUnset(pmStartSpan *sspan) { + if (!sspan) return false; + sspan->span = NULL; + sspan->direction = PM_STARTSPAN_NONE; + sspan->stop = false; + return true; +} + +static void pmFootprintSpansFree(pmFootprintSpans *fp) { + psFree(fp->startspans); + return; +} + +// Allocate a pmFootprintSpans structure with unassigned spans +pmFootprintSpans *pmFootprintSpansAlloc(int nSpans) { + pmFootprintSpans *fpSpans = psAlloc(sizeof(pmFootprintSpans)); + psMemSetDeallocator(fpSpans, (psFreeFunc)pmFootprintSpansFree); + + // create an footprint with allocated pmStartSpans, but none yet assigned + fpSpans->nStartSpans = 0; + fpSpans->startspans = psArrayAlloc(nSpans); + + for (int i = 0; i < nSpans; i++) { + fpSpans->startspans->data[i] = pmStartSpanAlloc(); + } + return fpSpans; +} + +// Allocate a pmFootprintSpans structure with unassigned spans +bool pmFootprintSpansInit(pmFootprintSpans *fpSpans) { + fpSpans->nStartSpans = 0; + for (int i = 0; i < fpSpans->startspans->n; i++) { + pmStartSpanUnset(fpSpans->startspans->data[i]); + } + return true; +} + +// Add a new pmSpan to a pmFootprintSpans. Iff we see a stop bit, return true +bool pmFootprintSpansSet(pmFootprintSpans *fpSpans, // the saved pmStartSpans + pmSpan *sp, // the span in question + psImage *mask, // mask of detected/stop pixels + const PM_STARTSPAN_DIR dir) { // the desired direction to search + if (dir == PM_STARTSPAN_RESTART) { + if (pmFootprintSpansSet(fpSpans, sp, mask, PM_STARTSPAN_UP)) return true; + if (pmFootprintSpansSet(fpSpans, sp, NULL, PM_STARTSPAN_DOWN)) return true; + } else { + int N = fpSpans->nStartSpans; + if (N == fpSpans->startspans->n) { + // if we need more space, extend fpSpans->startspans as needed + int Nalloc = fpSpans->startspans->n + 100; + psArrayRealloc(fpSpans->startspans, Nalloc); + fpSpans->startspans->n = Nalloc; + for (int i = N; i < fpSpans->startspans->n; i++) { + fpSpans->startspans->data[i] = pmStartSpanAlloc(); + } + } + pmStartSpan *startspan = fpSpans->startspans->data[N]; + + pmStartSpanSet (startspan, sp, mask, dir); + + if (startspan->stop) { // we detected a stop bit + pmStartSpanUnset(startspan); // don't allocate new span + return true; + } else { + fpSpans->nStartSpans ++; + return false; + } + } + return false; +} + Index: trunk/psModules/src/objects/pmFootprintSpans.h =================================================================== --- trunk/psModules/src/objects/pmFootprintSpans.h (revision 27672) +++ trunk/psModules/src/objects/pmFootprintSpans.h (revision 27672) @@ -0,0 +1,86 @@ +/* @file pmFootprintSpans.h + * + * @author RHL, Princeton & IfA; EAM, IfA + * + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-12-09 21:16:09 $ + * Copyright 2006 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_FOOTPRINT_SPANS_H +#define PM_FOOTPRINT_SPANS_H + +#include +#include "pmSpan.h" + +/* We define two helper structures used in building the pmFootprints: + * + * pmStartSpan : a smart span which knows how to (re-)start pixel scanning + * pmFootprintSpans : a collection of pmStartSpans which can define a footprint + * + * pmFootprintSpans allows us to allocate the memory for the spans before actually defining them + * + */ + +/* pmStartSpan + * + * A data structure to hold the starting point for a search for pixels above threshold, + * used by pmFootprintsFindAtPoint + * + * We don't want to find this span again --- it's already part of the footprint --- + * so we set appropriate mask bits + * + */ + +// +// An enum for what we should do with a pmStartSpan +// +typedef enum {PM_STARTSPAN_NONE = 0, // span is not defined + PM_STARTSPAN_DOWN, // scan down from this span + PM_STARTSPAN_UP, // scan up from this span + PM_STARTSPAN_RESTART, // restart scanning from this span + PM_STARTSPAN_DONE // this span is processed +} PM_STARTSPAN_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_STARTSPAN_STOP plane. It's simpler to be +// able to check for +// +enum { + PM_STARTSPAN_INITIAL = 0x0, // initial state of pixels. + PM_STARTSPAN_DETECTED = 0x1, // we've seen this pixel + PM_STARTSPAN_STOP = 0x2 // you may stop searching when you see this pixel +}; +// +// The struct that remembers how to [re-]start scanning the image for pixels +// +typedef struct { + pmSpan *span; // view on the real span (on a pmFootprint->spans array) + PM_STARTSPAN_DIR direction; // How to continue searching + bool stop; // should we stop searching? +} pmStartSpan; + +typedef struct { + psArray *startspans; + int nStartSpans; +} pmFootprintSpans; + +pmStartSpan *pmStartSpanAlloc(); + +bool pmStartSpanSet(pmStartSpan *sspan, + pmSpan *span, // The span in question + psImage *mask, // Pixels that we've already detected + const PM_STARTSPAN_DIR dir // Should we continue searching towards the top of the image? + ); + +pmFootprintSpans *pmFootprintSpansAlloc(int nSpans); + +bool pmFootprintSpansInit(pmFootprintSpans *fpSpans); + +bool pmFootprintSpansSet(pmFootprintSpans *fpSpans, // the saved pmStartSpans + pmSpan *sp, // the span in question + psImage *mask, // mask of detected/stop pixels + const PM_STARTSPAN_DIR dir); // the desired direction to search + +/// @} +# endif /* PM_FOOTPRINT_SPANS_H */ Index: trunk/psModules/src/psmodules.h =================================================================== --- trunk/psModules/src/psmodules.h (revision 27657) +++ trunk/psModules/src/psmodules.h (revision 27672) @@ -111,4 +111,5 @@ // the following headers are from psModule:objects #include +#include #include #include