Index: /trunk/psphot/src/pmFootprint.c =================================================================== --- /trunk/psphot/src/pmFootprint.c (revision 13441) +++ /trunk/psphot/src/pmFootprint.c (revision 13442) @@ -432,7 +432,9 @@ * * 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 - */ + * so we set appropriate mask bits + */ +// +// 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 @@ -440,352 +442,9 @@ 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 -} Old_Spartspan; - -static void Old_startspanFree(Old_Spartspan *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 Old_Spartspan * -Old_SpartspanAlloc(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? - ) { - Old_Spartspan *sspan = psAlloc(sizeof(Old_Spartspan)); - psMemSetDeallocator(sspan, (psFreeFunc)Old_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 Old_Spartspan to an array of Old_Spartspans -// -static void old_add_startspan(psArray *startspans, // the saved Old_Spartspans - 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) { - old_add_startspan(startspans, sp, img, PM_SSPAN_UP); - old_add_startspan(startspans, sp, NULL, PM_SSPAN_DOWN); - } else { - Old_Spartspan *sspan = Old_SpartspanAlloc(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 Old_Spartspan. - * 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 pmOldFindFootprintAtPoint - */ -static bool old_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; - - /********************************************************************************************************/ - - Old_Spartspan *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 Old_Spartspans 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); - - old_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); - old_add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_DONE); - } else { // overhangs to right - sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, x1 + col0); - old_add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_DONE); - sp = pmFootprintAddSpan(fp, i + row0, x1 + 1 + col0, sx1 + col0 - 1); - old_add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART); - } - first = false; - } else { - sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1); - old_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 * -pmOldFindFootprintAtPoint(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); - - old_add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART); - } - /* - * Now workout from those Old_Spartspans, searching for pixels above threshold - */ - while (old_do_startspan(fp, img, threshold, startspans)) continue; - /* - * Cleanup - */ - psFree(startspans); // restores the image pixel - - return fp; // pmFootprint really -} - -/************************************************************************************************************/ -/* - * 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 - */ -#if 0 -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 -#endif - // // 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. @@ -793,5 +452,7 @@ 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 @@ -1440,102 +1101,4 @@ * starting point, discard the peak. */ -psErrorCode pmOldFootprintCullPeaks(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 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); - - 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) - 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); - // - // 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; - 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 - (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]) { - threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON; - } - -#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 = pmOldFindFootprintAtPoint(subImg, threshold, peak->y, peak->x); - psImage *idImg = pmSetFootprintID(peakFootprint, peak_id); - psFree(peakFootprint); -#endif - - 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); - } - - psFree((psImage *)subImg); - psFree((psImage *)subWt); - - return PS_ERR_NONE; -} - -/************************************************************************************************************/ - /* - * 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(const psImage *img, // the image wherein lives the footprint const psImage *weight, // corresponding variance image @@ -1561,5 +1124,8 @@ // // We need a psArray of peaks brighter than the current peak. We'll fake this - // for efficiency --- shhh + // 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); @@ -1595,5 +1161,5 @@ } - const int peak_id = 1; // the ID for the peak of interest + const int peak_id = 1; // the ID for the peak of interest brightPeaks->n = i; // only stop at a peak brighter than we are pmFootprint *peakFootprint = pmFindFootprintAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x); Index: /trunk/psphot/src/pmFootprint.h =================================================================== --- /trunk/psphot/src/pmFootprint.h (revision 13441) +++ /trunk/psphot/src/pmFootprint.h (revision 13442) @@ -31,7 +31,4 @@ void pmFootprintSetBBox(pmFootprint *fp); psArray *pmFindFootprints(const psImage *img, const float threshold, const int npixMin); -pmFootprint *pmOldFindFootprintAtPoint(const psImage *img, - const float threshold, - int row, int col); pmFootprint *pmFindFootprintAtPoint(const psImage *img, const float threshold,