Index: trunk/psphot/src/Makefile.am =================================================================== --- trunk/psphot/src/Makefile.am (revision 17396) +++ trunk/psphot/src/Makefile.am (revision 17443) @@ -20,4 +20,13 @@ psphotErrorCodes.c \ pmFootprint.c \ + pmFootprintArrayGrow.c \ + pmFootprintArraysMerge.c \ + pmFootprintIDs.c \ + pmFootprintsFind.c \ + pmFootprintsFindAtPoint.c \ + pmFootprintsAssignPeaks.c \ + pmDetections.c \ + pmSpan.c \ + psphotCullPeaks.c \ psphotVersion.c \ psphotModelGroupInit.c \ @@ -29,4 +38,6 @@ psphotFindDetections.c \ psphotFindPeaks.c \ + psphotFindFootprints.c \ + psphotSignificanceImage.c \ psphotSourceStats.c \ psphotRoughClass.c \ Index: trunk/psphot/src/pmDetections.c =================================================================== --- trunk/psphot/src/pmDetections.c (revision 17443) +++ trunk/psphot/src/pmDetections.c (revision 17443) @@ -0,0 +1,26 @@ +# include "psphotInternal.h" + +void pmDetectionsFree (pmDetections *detections) { + + if (!detections) return; + + psFree (detections->footprints); + psFree (detections->peaks); + psFree (detections->oldPeaks); + return; +} + +// generate a pmDetections container with empty (allocated) footprints and peaks containers +pmDetections *pmDetectionsAlloc() { + + pmDetections *detections = (pmDetections *)psAlloc(sizeof(pmDetections)); + psMemSetDeallocator(detections, (psFreeFunc) pmDetectionsFree); + + detections->footprints = NULL; + detections->peaks = NULL; + detections->oldPeaks = NULL; + detections->last = 0; + + return (detections); +} + Index: trunk/psphot/src/pmFootprint.c =================================================================== --- trunk/psphot/src/pmFootprint.c (revision 17396) +++ trunk/psphot/src/pmFootprint.c (revision 17443) @@ -1,66 +1,5 @@ -/*****************************************************************************/ -/* - * Handle pmFootprints - */ -#include -#include -#include -#include -#include "psphot.h" -#include "pmFootprint.h" +# include "psphotInternal.h" -static void spanFree(pmSpan *tmp) { - ; -} - -/* - * pmSpanAlloc() - */ -pmSpan *pmSpanAlloc(int y, int x0, int x1) -{ - pmSpan *span = (pmSpan *)psAlloc(sizeof(pmSpan)); - psMemSetDeallocator(span, (psFreeFunc) spanFree); - - span->y = y; - span->x0 = x0; - span->x1 = x1; - - return(span); -} - -bool pmSpanTest(const psPtr ptr) -{ - return (psMemGetDeallocator(ptr) == (psFreeFunc)spanFree); -} - -// -// 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; - } -} - -/************************************************************************************************************/ +/**** low-level pmFootprint functions ****/ static void footprintFree(pmFootprint *tmp) @@ -131,8 +70,8 @@ // Add a span to a footprint, returning the new span // -static pmSpan *pmFootprintAddSpan(pmFootprint *fp, // the footprint to add to - const int y, // row to add - int x0, // range of - int x1) { // columns +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) { @@ -201,1020 +140,4 @@ } -/************************************************************************************************************/ - -typedef struct { /* run-length code for part of object*/ - int id; /* ID for object */ - int y; /* Row wherein WSPAN dwells */ - int x0, x1; /* inclusive range of columns */ -} WSPAN; - -/* - * comparison function for qsort; sort by ID then row - */ -static int -compar(const void *va, const void *vb) -{ - const WSPAN *a = va; - const WSPAN *b = vb; - - if(a->id < b->id) { - return(-1); - } else if(a->id > b->id) { - return(1); - } else { - return(a->y - b->y); - } -} - -/* - * Follow a chain of aliases, returning the final resolved value. - */ -static int -resolve_alias(const int *aliases, /* list of aliases */ - int id) /* alias to look up */ -{ - int resolved = id; /* resolved alias */ - - while(id != aliases[id]) { - resolved = id = aliases[id]; - } - - return(resolved); -} - -/* - * Go through an image, finding sets of connected pixels above threshold - * and assembling them into pmFootprints; the resulting set of objects - * is returned as a psArray - */ -psArray * -pmFindFootprints(const psImage *img, // image to search - const float threshold, // Threshold - const int npixMin) // minimum number of pixels in an acceptable pmFootprint -{ - int i0; /* initial value of i */ - int id; /* object ID */ - int in_span; /* object ID of current WSPAN */ - int nspan = 0; /* number of spans */ - int nobj = 0; /* number of objects found */ - int x0 = 0; /* unpacked from a WSPAN */ - int *tmp; /* used in swapping idc/idp */ - - 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; -/* - * Storage for arrays that identify objects by ID. We want to be able to - * refer to idp[-1] and idp[numCols], hence the (numCols + 2) - */ - int *id_s = psAlloc(2*(numCols + 2)*sizeof(int)); - memset(id_s, '\0', 2*(numCols + 2)*sizeof(int)); assert(id_s[0] == 0); - 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 - */ - for (int i = 0; i < numRows; i++) { - int j; - tmp = idc; idc = idp; idp = tmp; /* swap ID pointers */ - memset(idc, '\0', numCols*sizeof(int)); - - imgRowF32 = img->data.F32[i]; // only one of - imgRowS32 = img->data.S32[i]; // these is valid! - - in_span = 0; /* not in a span */ - for (j = 0; j < numCols; j++) { - double pixVal = F32 ? imgRowF32[j] : imgRowS32[j]; - if (pixVal < threshold) { - if (in_span) { - if(nspan >= size_spans) { - size_spans *= 2; - spans = psRealloc(spans, size_spans*sizeof(WSPAN)); - } - spans[nspan].id = in_span; - spans[nspan].y = i; - spans[nspan].x0 = x0; - spans[nspan].x1 = j - 1; - - nspan++; - - in_span = 0; - } - } else { /* a pixel to fix */ - if(idc[j - 1] != 0) { - id = idc[j - 1]; - } else if(idp[j - 1] != 0) { - id = idp[j - 1]; - } else if(idp[j] != 0) { - id = idp[j]; - } else if(idp[j + 1] != 0) { - id = idp[j + 1]; - } else { - id = ++nobj; - - if(id >= size_aliases) { - size_aliases *= 2; - aliases = psRealloc(aliases, size_aliases*sizeof(int)); - } - aliases[id] = id; - } - - idc[j] = id; - if(!in_span) { - x0 = j; in_span = id; - } -/* - * Do we need to merge ID numbers? If so, make suitable entries in aliases[] - */ - if(idp[j + 1] != 0 && idp[j + 1] != id) { - aliases[resolve_alias(aliases, idp[j + 1])] = - resolve_alias(aliases, id); - - idc[j] = id = idp[j + 1]; - } - } - } - - if(in_span) { - if(nspan >= size_spans) { - size_spans *= 2; - spans = psRealloc(spans, size_spans*sizeof(WSPAN)); - } - - assert(nspan < size_spans); /* we checked for space above */ - spans[nspan].id = in_span; - spans[nspan].y = i; - spans[nspan].x0 = x0; - spans[nspan].x1 = j - 1; - - nspan++; - } - } - - psFree(id_s); -/* - * Resolve aliases; first alias chains, then the IDs in the spans - */ - for (int i = 0; i < nspan; i++) { - spans[i].id = resolve_alias(aliases, spans[i].id); - } - - psFree(aliases); -/* - * Sort spans by ID, so we can sweep through them once - */ - if(nspan > 0) { - qsort(spans, nspan, sizeof(WSPAN), compar); - } -/* - * Build pmFootprints from the spans - */ - psArray *footprints = psArrayAlloc(nobj); - int n = 0; // number of pmFootprints - - if(nspan > 0) { - id = spans[0].id; - i0 = 0; - for (int i = 0; i <= nspan; i++) { /* nspan + 1 to catch last object */ - if(i == nspan || spans[i].id != id) { - pmFootprint *fp = pmFootprintAlloc(i - i0, img); - - for(; i0 < i; i0++) { - pmFootprintAddSpan(fp, spans[i0].y + row0, - spans[i0].x0 + col0, spans[i0].x1 + col0); - } - - if (fp->npix < npixMin) { - psFree(fp); - } else { - footprints->data[n++] = fp; - } - } - - id = spans[i].id; - } - } - - footprints = psArrayRealloc(footprints, n); -/* - * clean up - */ - psFree(spans); - - return footprints; -} - -/************************************************************************************************************/ -/* - * 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 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 - 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? -} Spartspan; - -static void startspanFree(Spartspan *sspan) { - psFree((void *)sspan->span); -} - -static Spartspan * -SpartspanAlloc(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? - ) { - Spartspan *sspan = psAlloc(sizeof(Spartspan)); - 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 - psMaskType *mpix = &mask->data.PS_TYPE_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 Spartspan to an array of Spartspans. Iff we see a stop bit, return true -// -static bool add_startspan(psArray *startspans, // the saved Spartspans - 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 { - Spartspan *sspan = SpartspanAlloc(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; -} - -/************************************************************************************************************/ -/* - * Search the image for pixels above threshold, starting at a single 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 pmFindFootprintAtPoint - */ -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? - 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 - psMaskType *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; - - /********************************************************************************************************/ - - Spartspan *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; - } - } - if (sspan == NULL || sspan->direction == PM_SSPAN_DONE) { // no more Spartspans to process - return false; - } - if (sspan->stop) { // they don't want any more spans processed - 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 - 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_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; - } - /* - * Cleanup - */ - - sspan->direction = PM_SSPAN_DONE; - return stop ? false : 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 - 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) { - 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_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_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! - psMaskType *maskRow = mask->data.PS_TYPE_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 Spartspans, 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 -} - -/************************************************************************************************************/ -/* - * 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 - int id = 0; // first index will be 1 - for (int i = 0; i < footprints->n; i++) { - const pmFootprint *fp = footprints->data[i]; - if (relativeIDs) { - id++; - } else { - id = fp->id; - } - - set_footprint_id(idImage, fp, id); - } -} - -/* - * Set an image to the value of footprint's ID whever they may fall - */ -psImage *pmSetFootprintArrayIDs(const psArray *footprints, // the footprints to insert - const bool relativeIDs) { // show IDs starting at 1, not pmFootprint->id - assert (footprints != NULL); - - if (footprints->n == 0) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, "You didn't provide any footprints"); - return NULL; - } - const pmFootprint *fp = footprints->data[0]; - assert(pmFootprintTest((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_array_ids(idImage, footprints, relativeIDs); - - return idImage; - -} - -/* - * 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 && pmFootprintTest((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; - -} - -/************************************************************************************************************/ -/* - * Grow a psArray of pmFootprints isotropically by r pixels, returning a new psArray of new pmFootprints - */ -psArray *pmGrowFootprintArray(const psArray *footprints, // footprints to grow - int r) { // how much to grow each footprint - assert (footprints->n == 0 || pmFootprintTest(footprints->data[0])); - - if (footprints->n == 0) { // we don't know the size of the footprint's region - return psArrayAlloc(0); - } - /* - * We'll insert the footprints into an image, then convolve with a disk, - * then extract a footprint from the result --- this is magically what we want. - */ - psImage *idImage = pmSetFootprintArrayIDs(footprints, true); - if (r <= 0) { - r = 1; // r == 1 => no grow - } - psKernel *circle = psKernelAlloc(-r, r, -r, r); - assert (circle->image->numRows == 2*r + 1 && circle->image->numCols == circle->image->numRows); - for (int i = 0; i <= r; i++) { - for (int j = 0; j <= r; j++) { - if (i*i + j*j <= r*r) { - circle->kernel[i][j] = - circle->kernel[i][-j] = - circle->kernel[-i][j] = - circle->kernel[-i][-j] = 1; - } - } - } - - psImage *grownIdImage = psImageConvolveDirect(NULL, idImage, circle); // Here's the actual grow step - psFree(circle); - - psArray *grown = pmFindFootprints(grownIdImage, 0.5, 1); // and here we rebuild the grown footprints - assert (grown != NULL); - psFree(idImage); psFree(grownIdImage); - /* - * Now assign the peaks appropriately. We could do this more efficiently - * using grownIdImage (which we just freed), but this is easy and probably fast enough - */ - const psArray *peaks = pmFootprintArrayToPeaks(footprints); - pmPeaksAssignToFootprints(grown, peaks); - psFree((psArray *)peaks); - - return grown; - -} - -/************************************************************************************************************/ -/* - * Merge together two psArrays of pmFootprints neither of which is damaged. - * - * The returned psArray may contain elements of the inital psArrays (with - * their reference counters suitable incremented) - */ -psArray *pmMergeFootprintArrays(const psArray *footprints1, // one set of footprints - const psArray *footprints2, // the other set - const int includePeaks) { // which peaks to set? 0x1 => footprints1, 0x2 => 2 - assert (footprints1->n == 0 || pmFootprintTest(footprints1->data[0])); - assert (footprints2->n == 0 || pmFootprintTest(footprints2->data[0])); - - if (footprints1->n == 0 || footprints2->n == 0) { // nothing to do but put copies on merged - const psArray *old = (footprints1->n == 0) ? footprints2 : footprints1; - - psArray *merged = psArrayAllocEmpty(old->n); - for (int i = 0; i < old->n; i++) { - psArrayAdd(merged, 1, old->data[i]); - } - - return merged; - } - /* - * We have real work to do as some pmFootprints in footprints2 may overlap - * with footprints1 - */ - { - pmFootprint *fp1 = footprints1->data[0]; - pmFootprint *fp2 = footprints2->data[0]; - if (fp1->region.x0 != fp2->region.x0 || - fp1->region.x1 != fp2->region.x1 || - fp1->region.y0 != fp2->region.y0 || - fp1->region.y1 != fp2->region.y1) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, - "The two pmFootprint arrays correspnond to different-sized regions"); - return NULL; - } - } - /* - * We'll insert first one set of footprints then the other into an image, then - * extract a footprint from the result --- this is magically what we want. - */ - psImage *idImage = pmSetFootprintArrayIDs(footprints1, true); - set_footprint_array_ids(idImage, footprints2, true); - - psArray *merged = pmFindFootprints(idImage, 0.5, 1); - assert (merged != NULL); - psFree(idImage); - /* - * Now assign the peaks appropriately. We could do this more efficiently - * using idImage (which we just freed), but this is easy and probably fast enough - */ - if (includePeaks & 0x1) { - const psArray *peaks = pmFootprintArrayToPeaks(footprints1); - pmPeaksAssignToFootprints(merged, peaks); - psFree((psArray *)peaks); - } - - if (includePeaks & 0x2) { - const psArray *peaks = pmFootprintArrayToPeaks(footprints2); - pmPeaksAssignToFootprints(merged, peaks); - psFree((psArray *)peaks); - } - - return merged; -} - -/************************************************************************************************************/ -/* - * Given a psArray of pmFootprints and another of pmPeaks, assign the peaks to the - * footprints in which that fall; if they _don't_ fall in a footprint, add a suitable - * one to the list. - */ -psErrorCode -pmPeaksAssignToFootprints(psArray *footprints, // the pmFootprints - const psArray *peaks) { // the pmPeaks - assert (footprints != NULL); - assert (footprints->n == 0 || pmFootprintTest(footprints->data[0])); - assert (peaks != NULL); - assert (peaks->n == 0 || psMemCheckPeak(peaks->data[0])); - - if (footprints->n == 0) { - if (peaks->n > 0) { - return psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Your list of footprints is empty"); - } - return PS_ERR_NONE; - } - /* - * Create an image filled with the object IDs, and use it to assign pmPeaks to the - * objects - */ - psImage *ids = pmSetFootprintArrayIDs(footprints, true); - assert (ids != NULL); - assert (ids->type.type == PS_TYPE_S32); - const int row0 = ids->row0; - const int col0 = ids->col0; - const int numRows = ids->numRows; - const int numCols = ids->numCols; - - for (int i = 0; i < peaks->n; i++) { - pmPeak *peak = peaks->data[i]; - const int x = peak->x - col0; - const int y = peak->y - row0; - - assert (x >= 0 && x < numCols && y >= 0 && y < numRows); - int id = ids->data.S32[y][x - col0]; - - if (id == 0) { // peak isn't in a footprint, so make one for it - pmFootprint *nfp = pmFootprintAlloc(1, ids); - pmFootprintAddSpan(nfp, y, x, x); - psArrayAdd(footprints, 1, nfp); - psFree(nfp); - id = footprints->n; - } - - assert (id >= 1 && id <= footprints->n); - pmFootprint *fp = footprints->data[id - 1]; - psArrayAdd(fp->peaks, 5, peak); - } - - psFree(ids); - // - // Make sure that peaks within each footprint are sorted and unique - // - for (int i = 0; i < footprints->n; i++) { - pmFootprint *fp = footprints->data[i]; - fp->peaks = psArraySort(fp->peaks, pmPeakSortBySN); - - for (int j = 1; j < fp->peaks->n; j++) { // check for duplicates - if (fp->peaks->data[j] == fp->peaks->data[j-1]) { - (void)psArrayRemoveIndex(fp->peaks, j); - j--; // we moved everything down one - } - } - } - - 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 - 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); - // - // 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); - 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; - } - - 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); - brightPeaks->n = 0; // don't double free - psImage *idImg = pmSetFootprintID(peakFootprint, peak_id); - psFree(peakFootprint); - - 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); - - return PS_ERR_NONE; -} - -/* - * Cull an entire psArray of pmFootprints - */ -psErrorCode -pmFootprintArrayCullPeaks(const psImage *img, // the image wherein lives the footprint - const psImage *weight, // corresponding variance image - psArray *footprints, // array of pmFootprints - 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 - for (int i = 0; i < footprints->n; i++) { - pmFootprint *fp = footprints->data[i]; - if (pmFootprintCullPeaks(img, weight, fp, nsigma_delta, min_threshold) != PS_ERR_NONE) { - return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id); - } - } - - return PS_ERR_NONE; -} - -/************************************************************************************************************/ /* * Extract the peaks in a psArray of pmFootprints, returning a psArray of pmPeaks @@ -1242,50 +165,3 @@ } -void pmDetectionsFree (pmDetections *detections) { - - if (!detections) return; - - psFree (detections->footprints); - psFree (detections->peaks); - psFree (detections->oldPeaks); - return; -} - -// generate a pmDetections container with empty (allocated) footprints and peaks containers -pmDetections *pmDetectionsAlloc() { - - pmDetections *detections = (pmDetections *)psAlloc(sizeof(pmDetections)); - psMemSetDeallocator(detections, (psFreeFunc) pmDetectionsFree); - - detections->footprints = NULL; - detections->peaks = NULL; - detections->oldPeaks = NULL; - detections->last = 0; - - return (detections); -} - /************************************************************************************************************/ -/* - * Cull a set of peaks contained in a psArray of pmFootprints - */ -psErrorCode -psphotCullPeaks(const psImage *image, // the image wherein lives the footprint - const psImage *weight, // corresponding variance image - const psMetadata *recipe, - psArray *footprints) { // array of pmFootprints - bool status = false; - float nsigma_delta = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_DELTA"); - if (!status) { - nsigma_delta = 0; // min. - } - float nsigma_min = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_MIN"); - if (!status) { - nsigma_min = 0; - } - const float skyStdev = psMetadataLookupF32(NULL, recipe, "SKY_STDEV"); - - return pmFootprintArrayCullPeaks(image, weight, footprints, - nsigma_delta, nsigma_min*skyStdev); -} - Index: trunk/psphot/src/pmFootprint.h =================================================================== --- trunk/psphot/src/pmFootprint.h (revision 17396) +++ trunk/psphot/src/pmFootprint.h (revision 17443) @@ -30,18 +30,27 @@ int pmFootprintSetNpix(pmFootprint *fp); void pmFootprintSetBBox(pmFootprint *fp); -psArray *pmFindFootprints(const psImage *img, const float threshold, const int npixMin); -pmFootprint *pmFindFootprintAtPoint(const psImage *img, + +pmSpan *pmFootprintAddSpan(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); -psArray *pmGrowFootprintArray(const psArray *footprints, int r); -psArray *pmMergeFootprintArrays(const psArray *footprints1, const psArray *footprints2, +psArray *pmFootprintArrayGrow(const psArray *footprints, int r); +psArray *pmFootprintArraysMerge(const psArray *footprints1, const psArray *footprints2, const int includePeaks); psImage *pmSetFootprintArrayIDs(const psArray *footprints, const bool relativeIDs); psImage *pmSetFootprintID(const pmFootprint *fp, const int id); +void pmSetFootprintArrayIDsForImage(psImage *idImage, + const psArray *footprints, // the footprints to insert + const bool relativeIDs); // show IDs starting at 0, not pmFootprint->id -psErrorCode pmPeaksAssignToFootprints(psArray *footprints, const psArray *peaks); +psErrorCode pmFootprintsAssignPeaks(psArray *footprints, const psArray *peaks); psErrorCode pmFootprintArrayCullPeaks(const psImage *img, const psImage *weight, psArray *footprints, Index: trunk/psphot/src/pmFootprintArrayGrow.c =================================================================== --- trunk/psphot/src/pmFootprintArrayGrow.c (revision 17443) +++ trunk/psphot/src/pmFootprintArrayGrow.c (revision 17443) @@ -0,0 +1,51 @@ +# include "psphotInternal.h" + +/* + * Grow a psArray of pmFootprints isotropically by r pixels, returning a new psArray of new pmFootprints + */ +psArray *pmFootprintArrayGrow(const psArray *footprints, // footprints to grow + int r) { // how much to grow each footprint + assert (footprints->n == 0 || pmFootprintTest(footprints->data[0])); + + if (footprints->n == 0) { // we don't know the size of the footprint's region + return psArrayAlloc(0); + } + /* + * We'll insert the footprints into an image, then convolve with a disk, + * then extract a footprint from the result --- this is magically what we want. + */ + psImage *idImage = pmSetFootprintArrayIDs(footprints, true); + if (r <= 0) { + r = 1; // r == 1 => no grow + } + psKernel *circle = psKernelAlloc(-r, r, -r, r); + assert (circle->image->numRows == 2*r + 1 && circle->image->numCols == circle->image->numRows); + for (int i = 0; i <= r; i++) { + for (int j = 0; j <= r; j++) { + if (i*i + j*j <= r*r) { + circle->kernel[i][j] = + circle->kernel[i][-j] = + circle->kernel[-i][j] = + circle->kernel[-i][-j] = 1; + } + } + } + + psImage *grownIdImage = psImageConvolveDirect(NULL, idImage, circle); // Here's the actual grow step + psFree(circle); + + psArray *grown = pmFootprintsFind(grownIdImage, 0.5, 1); // and here we rebuild the grown footprints + assert (grown != NULL); + psFree(idImage); psFree(grownIdImage); + /* + * Now assign the peaks appropriately. We could do this more efficiently + * using grownIdImage (which we just freed), but this is easy and probably fast enough + */ + const psArray *peaks = pmFootprintArrayToPeaks(footprints); + pmFootprintsAssignPeaks(grown, peaks); + psFree((psArray *)peaks); + + return grown; + +} + Index: trunk/psphot/src/pmFootprintArraysMerge.c =================================================================== --- trunk/psphot/src/pmFootprintArraysMerge.c (revision 17443) +++ trunk/psphot/src/pmFootprintArraysMerge.c (revision 17443) @@ -0,0 +1,68 @@ +# include "psphotInternal.h" + +/* + * Merge together two psArrays of pmFootprints neither of which is damaged. + * + * The returned psArray may contain elements of the inital psArrays (with + * their reference counters suitable incremented) + */ +psArray *pmFootprintArraysMerge(const psArray *footprints1, // one set of footprints + const psArray *footprints2, // the other set + const int includePeaks) { // which peaks to set? 0x1 => footprints1, 0x2 => 2 + assert (footprints1->n == 0 || pmFootprintTest(footprints1->data[0])); + assert (footprints2->n == 0 || pmFootprintTest(footprints2->data[0])); + + if (footprints1->n == 0 || footprints2->n == 0) { // nothing to do but put copies on merged + const psArray *old = (footprints1->n == 0) ? footprints2 : footprints1; + + psArray *merged = psArrayAllocEmpty(old->n); + for (int i = 0; i < old->n; i++) { + psArrayAdd(merged, 1, old->data[i]); + } + + return merged; + } + /* + * We have real work to do as some pmFootprints in footprints2 may overlap + * with footprints1 + */ + { + pmFootprint *fp1 = footprints1->data[0]; + pmFootprint *fp2 = footprints2->data[0]; + if (fp1->region.x0 != fp2->region.x0 || + fp1->region.x1 != fp2->region.x1 || + fp1->region.y0 != fp2->region.y0 || + fp1->region.y1 != fp2->region.y1) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "The two pmFootprint arrays correspnond to different-sized regions"); + return NULL; + } + } + /* + * We'll insert first one set of footprints then the other into an image, then + * extract a footprint from the result --- this is magically what we want. + */ + psImage *idImage = pmSetFootprintArrayIDs(footprints1, true); + pmSetFootprintArrayIDsForImage(idImage, footprints2, true); + + psArray *merged = pmFootprintsFind(idImage, 0.5, 1); + assert (merged != NULL); + psFree(idImage); + /* + * Now assign the peaks appropriately. We could do this more efficiently + * using idImage (which we just freed), but this is easy and probably fast enough + */ + if (includePeaks & 0x1) { + const psArray *peaks = pmFootprintArrayToPeaks(footprints1); + pmFootprintsAssignPeaks(merged, peaks); + psFree((psArray *)peaks); + } + + if (includePeaks & 0x2) { + const psArray *peaks = pmFootprintArrayToPeaks(footprints2); + pmFootprintsAssignPeaks(merged, peaks); + psFree((psArray *)peaks); + } + + return merged; +} Index: trunk/psphot/src/pmFootprintIDs.c =================================================================== --- trunk/psphot/src/pmFootprintIDs.c (revision 17443) +++ trunk/psphot/src/pmFootprintIDs.c (revision 17443) @@ -0,0 +1,97 @@ +# include "psphotInternal.h" + +/**** functions to manipulate footprint IDs ****/ + + +/* + * 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; + } + } +} + +void pmSetFootprintArrayIDsForImage(psImage *idImage, + const psArray *footprints, // the footprints to insert + const bool relativeIDs) { // show IDs starting at 0, not pmFootprint->id + int id = 0; // first index will be 1 + for (int i = 0; i < footprints->n; i++) { + const pmFootprint *fp = footprints->data[i]; + if (relativeIDs) { + id++; + } else { + id = fp->id; + } + + set_footprint_id(idImage, fp, id); + } +} + +/* + * Set an image to the value of footprint's ID whever they may fall + */ +psImage *pmSetFootprintArrayIDs(const psArray *footprints, // the footprints to insert + const bool relativeIDs) { // show IDs starting at 1, not pmFootprint->id + assert (footprints != NULL); + + if (footprints->n == 0) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, "You didn't provide any footprints"); + return NULL; + } + const pmFootprint *fp = footprints->data[0]; + assert(pmFootprintTest((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 + */ + pmSetFootprintArrayIDsForImage(idImage, footprints, relativeIDs); + + return idImage; + +} + +/* + * 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 && pmFootprintTest((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; + +} + Index: trunk/psphot/src/pmFootprintsAssignPeaks.c =================================================================== --- trunk/psphot/src/pmFootprintsAssignPeaks.c (revision 17443) +++ trunk/psphot/src/pmFootprintsAssignPeaks.c (revision 17443) @@ -0,0 +1,73 @@ +# include "psphotInternal.h" + +/* + * Given a psArray of pmFootprints and another of pmPeaks, assign the peaks to the + * footprints in which that fall; if they _don't_ fall in a footprint, add a suitable + * one to the list. + */ +psErrorCode +pmFootprintsAssignPeaks(psArray *footprints, // the pmFootprints + const psArray *peaks) { // the pmPeaks + assert (footprints != NULL); + assert (footprints->n == 0 || pmFootprintTest(footprints->data[0])); + assert (peaks != NULL); + assert (peaks->n == 0 || psMemCheckPeak(peaks->data[0])); + + if (footprints->n == 0) { + if (peaks->n > 0) { + return psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Your list of footprints is empty"); + } + return PS_ERR_NONE; + } + /* + * Create an image filled with the object IDs, and use it to assign pmPeaks to the + * objects + */ + psImage *ids = pmSetFootprintArrayIDs(footprints, true); + assert (ids != NULL); + assert (ids->type.type == PS_TYPE_S32); + const int row0 = ids->row0; + const int col0 = ids->col0; + const int numRows = ids->numRows; + const int numCols = ids->numCols; + + for (int i = 0; i < peaks->n; i++) { + pmPeak *peak = peaks->data[i]; + const int x = peak->x - col0; + const int y = peak->y - row0; + + assert (x >= 0 && x < numCols && y >= 0 && y < numRows); + int id = ids->data.S32[y][x - col0]; + + if (id == 0) { // peak isn't in a footprint, so make one for it + pmFootprint *nfp = pmFootprintAlloc(1, ids); + pmFootprintAddSpan(nfp, y, x, x); + psArrayAdd(footprints, 1, nfp); + psFree(nfp); + id = footprints->n; + } + + assert (id >= 1 && id <= footprints->n); + pmFootprint *fp = footprints->data[id - 1]; + psArrayAdd(fp->peaks, 5, peak); + } + + psFree(ids); + // + // Make sure that peaks within each footprint are sorted and unique + // + for (int i = 0; i < footprints->n; i++) { + pmFootprint *fp = footprints->data[i]; + fp->peaks = psArraySort(fp->peaks, pmPeakSortBySN); + + for (int j = 1; j < fp->peaks->n; j++) { // check for duplicates + if (fp->peaks->data[j] == fp->peaks->data[j-1]) { + (void)psArrayRemoveIndex(fp->peaks, j); + j--; // we moved everything down one + } + } + } + + return PS_ERR_NONE; +} + Index: trunk/psphot/src/pmFootprintsFind.c =================================================================== --- trunk/psphot/src/pmFootprintsFind.c (revision 17443) +++ trunk/psphot/src/pmFootprintsFind.c (revision 17443) @@ -0,0 +1,227 @@ +# include "psphotInternal.h" + +/**** find footprints in an image (fast on large scale images) ****/ + +// XXX EAM : why use WSPAN in here rather than pmSpan? +typedef struct { /* run-length code for part of object*/ + int id; /* ID for object */ + int y; /* Row wherein WSPAN dwells */ + int x0, x1; /* inclusive range of columns */ +} WSPAN; + +/* + * comparison function for qsort; sort by ID then row + */ +static int +compar(const void *va, const void *vb) +{ + const WSPAN *a = va; + const WSPAN *b = vb; + + if(a->id < b->id) { + return(-1); + } else if(a->id > b->id) { + return(1); + } else { + return(a->y - b->y); + } +} + +/* + * Follow a chain of aliases, returning the final resolved value. + */ +static int +resolve_alias(const int *aliases, /* list of aliases */ + int id) /* alias to look up */ +{ + int resolved = id; /* resolved alias */ + + while(id != aliases[id]) { + resolved = id = aliases[id]; + } + + return(resolved); +} + +/* + * Go through an image, finding sets of connected pixels above threshold + * and assembling them into pmFootprints; the resulting set of objects + * is returned as a psArray + */ +psArray * +pmFootprintsFind(const psImage *img, // image to search + const float threshold, // Threshold + const int npixMin) // minimum number of pixels in an acceptable pmFootprint +{ + int i0; /* initial value of i */ + int id; /* object ID */ + int in_span; /* object ID of current WSPAN */ + int nspan = 0; /* number of spans */ + int nobj = 0; /* number of objects found */ + int x0 = 0; /* unpacked from a WSPAN */ + int *tmp; /* used in swapping idc/idp */ + + 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; +/* + * Storage for arrays that identify objects by ID. We want to be able to + * refer to idp[-1] and idp[numCols], hence the (numCols + 2) + */ + int *id_s = psAlloc(2*(numCols + 2)*sizeof(int)); + memset(id_s, '\0', 2*(numCols + 2)*sizeof(int)); assert(id_s[0] == 0); + 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 + */ + for (int i = 0; i < numRows; i++) { + int j; + tmp = idc; idc = idp; idp = tmp; /* swap ID pointers */ + memset(idc, '\0', numCols*sizeof(int)); + + imgRowF32 = img->data.F32[i]; // only one of + imgRowS32 = img->data.S32[i]; // these is valid! + + in_span = 0; /* not in a span */ + for (j = 0; j < numCols; j++) { + double pixVal = F32 ? imgRowF32[j] : imgRowS32[j]; + if (pixVal < threshold) { + if (in_span) { + if(nspan >= size_spans) { + size_spans *= 2; + spans = psRealloc(spans, size_spans*sizeof(WSPAN)); + } + spans[nspan].id = in_span; + spans[nspan].y = i; + spans[nspan].x0 = x0; + spans[nspan].x1 = j - 1; + + nspan++; + + in_span = 0; + } + } else { /* a pixel to fix */ + if(idc[j - 1] != 0) { + id = idc[j - 1]; + } else if(idp[j - 1] != 0) { + id = idp[j - 1]; + } else if(idp[j] != 0) { + id = idp[j]; + } else if(idp[j + 1] != 0) { + id = idp[j + 1]; + } else { + id = ++nobj; + + if(id >= size_aliases) { + size_aliases *= 2; + aliases = psRealloc(aliases, size_aliases*sizeof(int)); + } + aliases[id] = id; + } + + idc[j] = id; + if(!in_span) { + x0 = j; in_span = id; + } +/* + * Do we need to merge ID numbers? If so, make suitable entries in aliases[] + */ + if(idp[j + 1] != 0 && idp[j + 1] != id) { + aliases[resolve_alias(aliases, idp[j + 1])] = + resolve_alias(aliases, id); + + idc[j] = id = idp[j + 1]; + } + } + } + + if(in_span) { + if(nspan >= size_spans) { + size_spans *= 2; + spans = psRealloc(spans, size_spans*sizeof(WSPAN)); + } + + assert(nspan < size_spans); /* we checked for space above */ + spans[nspan].id = in_span; + spans[nspan].y = i; + spans[nspan].x0 = x0; + spans[nspan].x1 = j - 1; + + nspan++; + } + } + + psFree(id_s); +/* + * Resolve aliases; first alias chains, then the IDs in the spans + */ + for (int i = 0; i < nspan; i++) { + spans[i].id = resolve_alias(aliases, spans[i].id); + } + + psFree(aliases); +/* + * Sort spans by ID, so we can sweep through them once + */ + if(nspan > 0) { + qsort(spans, nspan, sizeof(WSPAN), compar); + } +/* + * Build pmFootprints from the spans + */ + psArray *footprints = psArrayAlloc(nobj); + int n = 0; // number of pmFootprints + + if(nspan > 0) { + id = spans[0].id; + i0 = 0; + for (int i = 0; i <= nspan; i++) { /* nspan + 1 to catch last object */ + if(i == nspan || spans[i].id != id) { + pmFootprint *fp = pmFootprintAlloc(i - i0, img); + + for(; i0 < i; i0++) { + pmFootprintAddSpan(fp, spans[i0].y + row0, + spans[i0].x0 + col0, spans[i0].x1 + col0); + } + + if (fp->npix < npixMin) { + psFree(fp); + } else { + footprints->data[n++] = fp; + } + } + + id = spans[i].id; + } + } + + footprints = psArrayRealloc(footprints, n); +/* + * clean up + */ + psFree(spans); + + return footprints; +} Index: trunk/psphot/src/pmFootprintsFindAtPoint.c =================================================================== --- trunk/psphot/src/pmFootprintsFindAtPoint.c (revision 17443) +++ trunk/psphot/src/pmFootprintsFindAtPoint.c (revision 17443) @@ -0,0 +1,390 @@ +# include "psphotInternal.h" + +/**** find footprints in a small image relative to a reference point ****/ + +/* + * 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((void *)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 + psMaskType *mpix = &mask->data.PS_TYPE_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; +} + +/************************************************************************************************************/ +/* + * 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 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 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 + psMaskType *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; + } + } + 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; + } + /* + * 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 + 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_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; + } + /* + * Cleanup + */ + + sspan->direction = PM_SSPAN_DONE; + return stop ? false : 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 pmFootprintsFind 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 * +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) { + 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_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_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! + psMaskType *maskRow = mask->data.PS_TYPE_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 +} Index: trunk/psphot/src/pmSpan.c =================================================================== --- trunk/psphot/src/pmSpan.c (revision 17443) +++ trunk/psphot/src/pmSpan.c (revision 17443) @@ -0,0 +1,55 @@ +# include "psphotInternal.h" + +/*** functions to manipulate image spans ***/ + +static void spanFree(pmSpan *tmp) { + return; +} + +/* + * pmSpanAlloc() + */ +pmSpan *pmSpanAlloc(int y, int x0, int x1) +{ + pmSpan *span = (pmSpan *)psAlloc(sizeof(pmSpan)); + psMemSetDeallocator(span, (psFreeFunc) spanFree); + + span->y = y; + span->x0 = x0; + span->x1 = x1; + + return(span); +} + +bool pmSpanTest(const psPtr ptr) +{ + return (psMemGetDeallocator(ptr) == (psFreeFunc)spanFree); +} + +// +// 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; + } +} Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 17396) +++ trunk/psphot/src/psphot.h (revision 17443) @@ -49,5 +49,7 @@ // used by psphotFindDetections -psArray *psphotFindPeaks (pmReadout *readout, psMetadata *recipe, const bool returnFootprints, const int pass, psMaskType maskVal); +psImage *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, const int pass, psMaskType maskVal); +psArray *psphotFindPeaks (psImage *significance, pmReadout *readout, psMetadata *recipe, const int pass, psMaskType maskVal); +bool psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const int pass, psMaskType maskVal); psErrorCode psphotCullPeaks(const psImage *img, const psImage *weight, const psMetadata *recipe, psArray *footprints); Index: trunk/psphot/src/psphotCullPeaks.c =================================================================== --- trunk/psphot/src/psphotCullPeaks.c (revision 17443) +++ trunk/psphot/src/psphotCullPeaks.c (revision 17443) @@ -0,0 +1,159 @@ +# include "psphotInternal.h" + +/* + * Cull a set of peaks contained in a psArray of pmFootprints + */ +psErrorCode +psphotCullPeaks(const psImage *image, // the image wherein lives the footprint + const psImage *weight, // corresponding variance image + const psMetadata *recipe, + psArray *footprints) { // array of pmFootprints + bool status = false; + float nsigma_delta = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_DELTA"); + if (!status) { + nsigma_delta = 0; // min. + } + float nsigma_min = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_MIN"); + if (!status) { + nsigma_min = 0; + } + const float skyStdev = psMetadataLookupF32(NULL, recipe, "SKY_STDEV"); + + return pmFootprintArrayCullPeaks(image, weight, footprints, + nsigma_delta, nsigma_min*skyStdev); +} + + +/* + * Cull an entire psArray of pmFootprints + * XXX drop this intermediate level function? + */ +psErrorCode +pmFootprintArrayCullPeaks(const psImage *img, // the image wherein lives the footprint + const psImage *weight, // corresponding variance image + psArray *footprints, // array of pmFootprints + 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 + for (int i = 0; i < footprints->n; i++) { + pmFootprint *fp = footprints->data[i]; + if (pmFootprintCullPeaks(img, weight, fp, nsigma_delta, min_threshold) != PS_ERR_NONE) { + return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id); + } + } + + 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 + 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); + // + // 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); + 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 + // 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(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); + + return PS_ERR_NONE; +} + Index: trunk/psphot/src/psphotFindFootprints.c =================================================================== --- trunk/psphot/src/psphotFindFootprints.c (revision 17443) +++ trunk/psphot/src/psphotFindFootprints.c (revision 17443) @@ -0,0 +1,79 @@ +# include "psphotInternal.h" + +// N.b. We're not culling this list; call pmFootprintCullPeaks if you +// want to do that + +bool psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const int pass, psMaskType maskVal) { + + bool status; + + psTimerStart ("footprints"); + + int npixMin = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NPIXMIN"); + PS_ASSERT (status, false); + + float FOOTPRINT_NSIGMA_LIMIT; + if (pass == 1) { + FOOTPRINT_NSIGMA_LIMIT = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NSIGMA_LIMIT"); + } else { + FOOTPRINT_NSIGMA_LIMIT = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NSIGMA_LIMIT_2"); + } + PS_ASSERT (status, false); + + float effArea = psMetadataLookupF32 (&status, recipe, "EFFECTIVE_AREA"); + psAssert (status, "EFFECTIVE_AREA missing: call psphotFindPeaks first"); + + float threshold = PS_SQR(FOOTPRINT_NSIGMA_LIMIT)/effArea; + + int growRadius = 0; + if (pass == 1) { + growRadius = psMetadataLookupS32(&status, recipe, "FOOTPRINT_GROW_RADIUS"); + } else { + growRadius = psMetadataLookupS32(&status, recipe, "FOOTPRINT_GROW_RADIUS_2"); + } + PS_ASSERT (status, false); + + // XXX do we need to use the same threshold here as for peaks? does it make sense for + // these to be different? + + // find the raw footprints & assign the peaks to those footprints + psArray *footprints = pmFootprintsFind (significance, threshold, npixMin); + pmFootprintsAssignPeaks(footprints, detections->peaks); + // XXX handle the error conditions here + + // footprints now owns the peaks; after culling (below), we will rebuild + // the peaks array + psFree (detections->peaks); + + psLogMsg ("psphot", PS_LOG_DETAIL, "found %ld footprints: %f sec\n", footprints->n, psTimerMark ("footprints")); + + // optionally grow footprints isotropically by growRadius pixels + if (growRadius > 0) { + psArray *tmp = pmFootprintArrayGrow(footprints, growRadius); + psLogMsg ("psphot", PS_LOG_DETAIL, "grow footprint coverage by %d pixels, %ld footprints become %ld footprints: %f sec\n", growRadius, footprints->n, tmp->n, psTimerMark ("footprints")); + psFree(footprints); + footprints = tmp; + } + + if (pass == 2) { + // merge in old peaks; + const int includePeaks = 0x0 | 0x2; // i.e. just from newFootprints + + // XXX EAM : still not sure I understand this: are we double-couning or undercounting peaks + + psArray *mergedFootprints = pmFootprintArraysMerge(detections->footprints, footprints, includePeaks); + psLogMsg ("psphot", PS_LOG_DETAIL, "merged %ld new footprints with %ld existing ones: %f sec\n", footprints->n, detections->footprints->n, psTimerMark ("footprints")); + + psFree(footprints); + psFree(detections->footprints); + detections->footprints = mergedFootprints; + } else { + detections->footprints = footprints; + } + + psphotCullPeaks(readout->image, readout->weight, recipe, detections->footprints); + detections->peaks = pmFootprintArrayToPeaks(detections->footprints); + psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("footprints")); + + return detections; +} Index: trunk/psphot/src/psphotSignificanceImage.c =================================================================== --- trunk/psphot/src/psphotSignificanceImage.c (revision 17443) +++ trunk/psphot/src/psphotSignificanceImage.c (revision 17443) @@ -0,0 +1,76 @@ +# include "psphotInternal.h" + +// In this function, we smooth the image, then search for the peaks +// if FWMH_X,Y have been recorded, use them; otherwise use PEAKS_SMOOTH_SIGMA +psImage *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, const int pass, psMaskType maskVal) { + + float SIGMA_SMTH, NSIGMA_SMTH; + bool status = false; + + // smooth the image and weight map + psTimerStart ("smooth"); + + bool status_x, status_y; + float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X"); + float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y"); + if (status_x && status_y) { + // if we know the FHWM, use that to set the smoothing kernel (XXX allow an optional override?) + SIGMA_SMTH = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrt(2.0*log(2.0))); + NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA"); + } else { + // if we do not know the FWHM, use the guess smoothing kernel supplied. + // it is a configuration error if these are not supplied + SIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA"); + PS_ASSERT (status, NULL); + NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA"); + PS_ASSERT (status, NULL); + } + // record the actual smoothing sigma + psMetadataAddF32 (recipe, PS_LIST_TAIL, "SIGMA_SMOOTH", PS_META_REPLACE, "Smoothing sigma for detections", SIGMA_SMTH); + + // smooth the image, applying the mask as we go + psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32); + psImageSmoothMaskF32 (smooth_im, readout->mask, maskVal, SIGMA_SMTH, NSIGMA_SMTH); + psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("smooth")); + + // smooth the weight, applying the mask as we go + psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32); + psImageSmoothMaskF32 (smooth_wt, readout->mask, maskVal, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH); + psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("smooth")); + + psImage *mask = readout->mask; + + // optionally save example images under trace + // XXX change these to recipe value checks + if (psTraceGetLevel("psphot") > 5) { + char name[64]; + sprintf (name, "imsmooth.v%d.fits", pass); + psphotSaveImage (NULL, smooth_im, name); + sprintf (name, "wtsmooth.v%d.fits", pass); + psphotSaveImage (NULL, smooth_wt, name); + } + + // build the significance image on top of smooth_im + for (int j = 0; j < smooth_im->numRows; j++) { + for (int i = 0; i < smooth_im->numCols; i++) { + float value = smooth_im->data.F32[j][i]; + if (value < 0 || smooth_wt->data.F32[j][i] <= 0 || (mask->data.U8[j][i] & maskVal)) { + smooth_im->data.F32[j][i] = 0.0; + } else { + smooth_im->data.F32[j][i] = PS_SQR(value) / smooth_wt->data.F32[j][i]; + } + } + } + psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark ("smooth")); + + // optionally save example images under trace + // XXX change these to recipe value checks + if (psTraceGetLevel("psphot") > 5) { + char name[64]; + sprintf (name, "snsmooth.v%d.fits", pass); + psphotSaveImage (NULL, smooth_im, name); + } + + psFree (smooth_wt); + return smooth_im; +}