Index: trunk/psphot/src/pmFootprint.c =================================================================== --- trunk/psphot/src/pmFootprint.c (revision 13007) +++ trunk/psphot/src/pmFootprint.c (revision 13236) @@ -34,4 +34,32 @@ } +// +// Sort pmSpans by y, then x0, then x1 +// +int pmSpanSortByYX(const void **a, const void **b) { + const pmSpan *sa = *(const pmSpan **)a; + const pmSpan *sb = *(const pmSpan **)b; + + if (sa->y < sb->y) { + return -1; + } else if (sa->y == sb->y) { + if (sa->x0 < sb->x0) { + return -1; + } else if (sa->x0 == sb->x0) { + if (sa->x1 < sb->x1) { + return -1; + } else if (sa->x1 == sb->x1) { + return 0; + } else { + return 1; + } + } else { + return 1; + } + } else { + return 1; + } +} + /************************************************************************************************************/ @@ -63,7 +91,8 @@ psMemSetDeallocator(footprint, (psFreeFunc) footprintFree); + footprint->normalized = false; + assert(nspan >= 0); footprint->npix = 0; - footprint->nspan = nspan; footprint->spans = psArrayAllocEmpty(nspan); footprint->peaks = psArrayAlloc(0); @@ -90,11 +119,20 @@ } +pmFootprint *pmFootprintNormalize(pmFootprint *fp) { + if (fp != NULL && !fp->normalized) { + fp->peaks = psArraySort(fp->peaks, pmPeakSortBySN); + fp->normalized = true; + } + + return fp; +} + // -// Add a span to a footprint, returning the number of pixels in the footprint +// Add a span to a footprint, returning the new span // -static int pmFootprintAddSpan(pmFootprint *fp, // the footprint to add to - const int y, // row to add - int x0, // range of - int x1) { // columns +static pmSpan *pmFootprintAddSpan(pmFootprint *fp, // the footprint to add to + const int y, // row to add + int x0, // range of + int x1) { // columns if (x1 < x0) { @@ -110,5 +148,5 @@ fp->npix += x1 - x0 + 1; - if (fp->nspan == 1) { + if (fp->spans->n == 1) { fp->bbox.x0 = x0; fp->bbox.x1 = x1; @@ -122,10 +160,10 @@ } - return fp->npix; + return sp; } void pmFootprintSetBBox(pmFootprint *fp) { assert (fp != NULL); - if (fp->nspan == 0) { + if (fp->spans->n == 0) { return; } @@ -154,5 +192,5 @@ assert (fp != NULL); int npix = 0; - for (int i = 0; i < fp->nspan; i++) { + for (int i = 0; i < fp->spans->n; i++) { pmSpan *span = fp->spans->data[i]; npix += span->x1 - span->x0 + 1; @@ -215,7 +253,4 @@ const int npixMin) // minimum number of pixels in an acceptable pmFootprint { - int *aliases; /* aliases for object IDs */ - int *id_s; /* storage for id[cp] */ - int *idc, *idp; /* object IDs in current/previous row*/ int i0; /* initial value of i */ int id; /* object ID */ @@ -224,7 +259,4 @@ int nobj = 0; /* number of objects found */ int x0 = 0; /* unpacked from a WSPAN */ - int size_aliases = 0; /* size of aliases[] array */ - int size_spans = 0; /* size of spans[] array */ - WSPAN *spans; /* row:x0,x1 for objects */ int *tmp; /* used in swapping idc/idp */ @@ -251,13 +283,14 @@ * refer to idp[-1] and idp[numCols], hence the (numCols + 2) */ - id_s = psAlloc(2*(numCols + 2)*sizeof(int)); + int *id_s = psAlloc(2*(numCols + 2)*sizeof(int)); memset(id_s, '\0', 2*(numCols + 2)*sizeof(int)); assert(id_s[0] == 0); - idc = id_s + 1; idp = idc + (numCols + 2); - - size_aliases = 1 + numRows/20; - aliases = psAlloc(size_aliases*sizeof(int)); - - size_spans = 1 + numRows/20; - spans = psAlloc(size_spans*sizeof(WSPAN)); + int *idc = id_s + 1; // object IDs in current/ + int *idp = idc + (numCols + 2); // previous row + + int size_aliases = 1 + numRows/20; // size of aliases[] array + int *aliases = psAlloc(size_aliases*sizeof(int)); // aliases for object IDs + + int size_spans = 1 + numRows/20; // size of spans[] array + WSPAN *spans = psAlloc(size_spans*sizeof(WSPAN)); // row:x0,x1 for objects /* * Go through image identifying objects @@ -273,5 +306,5 @@ in_span = 0; /* not in a span */ for (j = 0; j < numCols; j++) { - double pixVal = F32 ? imgRowF32[j] : imgRowS32[j]; + double pixVal = F32 ? imgRowF32[j] : imgRowS32[j]; if (pixVal < threshold) { if (in_span) { @@ -395,30 +428,377 @@ /************************************************************************************************************/ /* - * Worker routine for the pmSetFootprintArrayIds (and pmMergeFootprintArrays) - */ + * A data structure to hold the starting point for a search for pixels above threshold, + * used by pmFindFootprintAtPoint + * + * We don't want to find this span again --- it's already part of the footprint --- + * so we set the pixels in the image to as small a value as we can. N.b. We may want + * to reconsider this when dealing with difference images + */ +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 struct { + const pmSpan *span; // save the pixel range + PM_SSPAN_DIR direction; // How to continue searching + // private + int nbyte; // number of bytes saved + void *data; // saved data + void *data_home; // where the data was saved from +} Startspan; + +static void startspanFree(Startspan *sspan) { + if (sspan->nbyte > 0) { + memcpy(sspan->data_home, sspan->data, sspan->nbyte); // restore image's pixels to their pristine state + } + + psFree(sspan->data); + psFree((void *)sspan->span); +} + +static Startspan * +StartspanAlloc(const pmSpan *span, // The span in question + psImage *img,// the data wherein lives the span + 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; + + if (img == NULL) { // nothing to save + sspan->nbyte = 0; + sspan->data_home = sspan->data = NULL; + + return sspan; + } + // + // Save the pixels, and set the image to a -ve large value + // + sspan->nbyte = (span->x1 - span->x0 + 1)*PSELEMTYPE_SIZEOF(img->type.type); + sspan->data = psAlloc(sspan->nbyte); + + switch (img->type.type) { + case PS_TYPE_F32: + sspan->data_home = &img->data.F32[span->y - img->row0][span->x0 - img->col0]; + memcpy(sspan->data, sspan->data_home, sspan->nbyte); + + for (int i = 0; i <= span->x1 - span->x0; i++) { + ((psF32 *)sspan->data_home)[i] = PS_MIN_F32; + } + break; + case PS_TYPE_S32: + sspan->data_home = &img->data.S32[span->y - img->row0][span->x0 - img->col0]; + memcpy(sspan->data, sspan->data_home, sspan->nbyte); + + for (int i = 0; i <= span->x1 - span->x0; i++) { + ((psS32 *)sspan->data_home)[i] = PS_MIN_S32; + } + break; + default: + psAbort("Unsupported image type %d\n", img->type.type); + } + + return sspan; +} + +// +// Add a new Startspan to an array of Startspans +// +static void add_startspan(psArray *startspans, // the saved Startspans + const pmSpan *sp, // the span in question + const psImage *img, // image to save/restore (or NULL) + const PM_SSPAN_DIR dir) { // the desired direction to search + if (dir == PM_SSPAN_RESTART) { + add_startspan(startspans, sp, img, PM_SSPAN_UP); + add_startspan(startspans, sp, NULL, PM_SSPAN_DOWN); + } else { + Startspan *sspan = StartspanAlloc(sp, (psImage *)img, dir); + psArrayAdd(startspans, 1, sspan); + psFree(sspan); // as it's now owned by startspans + } +} + +/************************************************************************************************************/ +/* + * Search the image for pixels above threshold, starting at a single Startspan. + * 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 + * the anticipated uses of this routine. + * + * This is the guts of pmFindFootprintAtPoint + */ +static bool do_startspan(pmFootprint *fp, // the footprint that we're building + const psImage *img, // the psImage we're working on + 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; + } 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; + + /********************************************************************************************************/ + + Startspan *sspan = NULL; + for (int i = 0; i < startspans->n; i++) { + sspan = startspans->data[i]; + if (sspan->direction != PM_SSPAN_DONE) { + break; + } + } + if (sspan == NULL || sspan->direction == PM_SSPAN_DONE) { // no more Startspans to process + return false; + } + /* + * Work + */ + const PM_SSPAN_DIR dir = sspan->direction; + /* + * Set initial span to the startspan + */ + int x0 = sspan->span->x0 - col0, x1 = sspan->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 + 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! + // + // 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 (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 (pixVal < threshold) { + nx1 = j - 1; + break; + } + } + + const pmSpan *sp = pmFootprintAddSpan(fp, i + row0, nx0 + col0, nx1 + col0); + + add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART); + } + // + // 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 (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 (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); + add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_DONE); + } else { // overhangs to right + sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, x1 + col0); + add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_DONE); + sp = pmFootprintAddSpan(fp, i + row0, x1 + 1 + col0, sx1 + col0 - 1); + add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART); + } + first = false; + } else { + sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1); + add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART); + } + } + } + + if (first == false) { // we're done + break; + } + + x0 = nx0; x1 = nx1; + } + /* + * Cleanup + */ + + sspan->direction = PM_SSPAN_DONE; + return true; +} + +/* + * Go through an image, starting at (row, col) and assembling all the pixels + * that are connected to that point (in a chess kings-move sort of way) into + * a pmFootprint. + * + * This is much slower than pmFindFootprints if you want to find lots of + * footprints, but if you only want a small region about a given point it + * can be much faster + * + * N.b. The returned pmFootprint is not in "normal form"; that is the pmSpans + * are not sorted by increasing y, x0, x1. If this matters to you, call + * pmFootprintNormalize() + */ +pmFootprint * +pmFindFootprintAtPoint(const psImage *img, // image to search + const float threshold, // Threshold + 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) { + 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); +/* + * 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! + { + int i; + for (i = col; i >= 0; i--) { + pixVal = F32 ? imgRowF32[i] : imgRowS32[i]; + if (pixVal < threshold) { + break; + } + } + int i0 = i; + for (i = col; i < numCols; i++) { + pixVal = F32 ? imgRowF32[i] : imgRowS32[i]; + if (pixVal < threshold) { + break; + } + } + int i1 = i; + const pmSpan *sp = pmFootprintAddSpan(fp, row + row0, i0 + col0 + 1, i1 + col0 - 1); + + add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART); + } + /* + * Now workout from those Startspans, searching for pixels above threshold + */ + while (do_startspan(fp, img, threshold, startspans)) continue; + /* + * Cleanup + */ + psFree(startspans); // restores the image pixel + + return fp; // pmFootprint really +} + +/************************************************************************************************************/ +/* + * Worker routine for the pmSetFootprintArrayIDs/pmSetFootprintID (and pmMergeFootprintArrays) + */ +static void +set_footprint_id(psImage *idImage, // the image to set + const pmFootprint *fp, // the footprint to insert + const int id) { // the desired ID + const int col0 = fp->region.x0; + const int row0 = fp->region.y0; + + for (int j = 0; j < fp->spans->n; j++) { + const pmSpan *span = fp->spans->data[j]; + psS32 *imgRow = idImage->data.S32[span->y - row0]; + for(int k = span->x0 - col0; k <= span->x1 - col0; k++) { + imgRow[k] += id; + } + } +} + static void set_footprint_array_ids(psImage *idImage, const psArray *footprints, // the footprints to insert const bool relativeIDs) { // show IDs starting at 0, not pmFootprint->id - const pmFootprint *fp = footprints->data[0]; - const int col0 = fp->region.x0; - const int row0 = fp->region.y0; - int id = 0; // first index will be 1 for (int i = 0; i < footprints->n; i++) { - fp = footprints->data[i]; - if (relativeIDs) { - id++; - } else { - id = fp->id; - } - - for (int j = 0; j < fp->nspan; j++) { - const pmSpan *span = fp->spans->data[j]; - psS32 *imgRow = idImage->data.S32[span->y - row0]; - for(int k = span->x0 - col0; k <= span->x1 - col0; k++) { - imgRow[k] += id; - } - } + const pmFootprint *fp = footprints->data[i]; + if (relativeIDs) { + id++; + } else { + id = fp->id; + } + + set_footprint_id(idImage, fp, id); } } @@ -428,5 +808,5 @@ */ psImage *pmSetFootprintArrayIDs(const psArray *footprints, // the footprints to insert - const bool relativeIDs) { // show IDs starting at 0, not pmFootprint->id + const bool relativeIDs) { // show IDs starting at 1, not pmFootprint->id assert (footprints != NULL); @@ -456,4 +836,29 @@ } +/* + * Set an image to the value of footprint's ID whever they may fall + */ +psImage *pmSetFootprintID(const pmFootprint *fp, // the footprint to insert + const int id) { // the desired ID + assert(fp != NULL && pmIsFootprint((const psPtr)fp)); + const int numCols = fp->region.x1 - fp->region.x0 + 1; + const int numRows = fp->region.y1 - fp->region.y0 + 1; + const int col0 = fp->region.x0; + const int row0 = fp->region.y0; + assert (numCols >= 0 && numRows >= 0); + + psImage *idImage = psImageAlloc(numCols, numRows, PS_TYPE_S32); + P_PSIMAGE_SET_ROW0(idImage, row0); + P_PSIMAGE_SET_COL0(idImage, col0); + psImageInit(idImage, 0); + /* + * do the work + */ + set_footprint_id(idImage, fp, id); + + return idImage; + +} + /************************************************************************************************************/ /* @@ -575,5 +980,5 @@ for (int i = 0; i < footprints->n; i++) { pmFootprint *fp = footprints->data[i]; - fp->peaks = psArraySort (fp->peaks, pmPeakSortBySN); + fp->peaks = psArraySort(fp->peaks, pmPeakSortBySN); } @@ -624,5 +1029,13 @@ float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev; if (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 + (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 } if (threshold > subImg->data.F32[y][x]) { @@ -630,9 +1043,17 @@ } - psArray *peakFootprints = pmFindFootprints(subImg, threshold, 1); +#if 0 + const in npixMin = 1; + psArray *peakFootprints = pmFindFootprints(subImg, threshold, npixMin); psImage *idImg = pmSetFootprintArrayIDs(peakFootprints, true); - + psFree(peakFootprints); const int peak_id = idImg->data.S32[y][x]; // the ID for the peak of interest assert (peak_id != 0); +#else + const int peak_id = 1; // the ID for the peak of interest + pmFootprint *peakFootprint = pmFindFootprintAtPoint(subImg, threshold, peak->y, peak->x); + psImage *idImg = pmSetFootprintID(peakFootprint, peak_id); + psFree(peakFootprint); +#endif int j; @@ -654,5 +1075,4 @@ } - psFree(peakFootprints); psFree(idImg); }