Changeset 13442 for trunk/psphot/src/pmFootprint.c
- Timestamp:
- May 21, 2007, 8:49:24 AM (19 years ago)
- File:
-
- 1 edited
-
trunk/psphot/src/pmFootprint.c (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psphot/src/pmFootprint.c
r13441 r13442 432 432 * 433 433 * We don't want to find this span again --- it's already part of the footprint --- 434 * so we set the pixels in the image to as small a value as we can. N.b. We may want 435 * to reconsider this when dealing with difference images 436 */ 434 * so we set appropriate mask bits 435 */ 436 // 437 // An enum for what we should do with a Startspan 438 // 437 439 typedef enum {PM_SSPAN_DOWN = 0, // scan down from this span 438 440 PM_SSPAN_UP, // scan up from this span … … 440 442 PM_SSPAN_DONE // this span is processed 441 443 } PM_SSPAN_DIR; // How to continue searching 442 443 typedef struct {444 const pmSpan *span; // save the pixel range445 PM_SSPAN_DIR direction; // How to continue searching446 // private447 int nbyte; // number of bytes saved448 void *data; // saved data449 void *data_home; // where the data was saved from450 } Old_Spartspan;451 452 static void Old_startspanFree(Old_Spartspan *sspan) {453 if (sspan->nbyte > 0) {454 memcpy(sspan->data_home, sspan->data, sspan->nbyte); // restore image's pixels to their pristine state455 }456 457 psFree(sspan->data);458 psFree((void *)sspan->span);459 }460 461 static Old_Spartspan *462 Old_SpartspanAlloc(const pmSpan *span, // The span in question463 psImage *img,// the data wherein lives the span464 const PM_SSPAN_DIR dir // Should we continue searching towards the top of the image?465 ) {466 Old_Spartspan *sspan = psAlloc(sizeof(Old_Spartspan));467 psMemSetDeallocator(sspan, (psFreeFunc)Old_startspanFree);468 469 sspan->span = psMemIncrRefCounter((void *)span);470 sspan->direction = dir;471 472 if (img == NULL) { // nothing to save473 sspan->nbyte = 0;474 sspan->data_home = sspan->data = NULL;475 476 return sspan;477 }478 //479 // Save the pixels, and set the image to a -ve large value480 //481 sspan->nbyte = (span->x1 - span->x0 + 1)*PSELEMTYPE_SIZEOF(img->type.type);482 sspan->data = psAlloc(sspan->nbyte);483 484 switch (img->type.type) {485 case PS_TYPE_F32:486 sspan->data_home = &img->data.F32[span->y - img->row0][span->x0 - img->col0];487 memcpy(sspan->data, sspan->data_home, sspan->nbyte);488 489 for (int i = 0; i <= span->x1 - span->x0; i++) {490 ((psF32 *)sspan->data_home)[i] = PS_MIN_F32;491 }492 break;493 case PS_TYPE_S32:494 sspan->data_home = &img->data.S32[span->y - img->row0][span->x0 - img->col0];495 memcpy(sspan->data, sspan->data_home, sspan->nbyte);496 497 for (int i = 0; i <= span->x1 - span->x0; i++) {498 ((psS32 *)sspan->data_home)[i] = PS_MIN_S32;499 }500 break;501 default:502 psAbort("Unsupported image type %d\n", img->type.type);503 }504 505 return sspan;506 }507 508 //509 // Add a new Old_Spartspan to an array of Old_Spartspans510 //511 static void old_add_startspan(psArray *startspans, // the saved Old_Spartspans512 const pmSpan *sp, // the span in question513 const psImage *img, // image to save/restore (or NULL)514 const PM_SSPAN_DIR dir) { // the desired direction to search515 if (dir == PM_SSPAN_RESTART) {516 old_add_startspan(startspans, sp, img, PM_SSPAN_UP);517 old_add_startspan(startspans, sp, NULL, PM_SSPAN_DOWN);518 } else {519 Old_Spartspan *sspan = Old_SpartspanAlloc(sp, (psImage *)img, dir);520 psArrayAdd(startspans, 1, sspan);521 psFree(sspan); // as it's now owned by startspans522 }523 }524 525 /************************************************************************************************************/526 /*527 * Search the image for pixels above threshold, starting at a single Old_Spartspan.528 * We search the array looking for one to process; it'd be better to move the529 * ones that we're done with to the end, but it probably isn't worth it for530 * the anticipated uses of this routine.531 *532 * This is the guts of pmOldFindFootprintAtPoint533 */534 static bool old_do_startspan(pmFootprint *fp, // the footprint that we're building535 const psImage *img, // the psImage we're working on536 const float threshold, // Threshold537 psArray *startspans) { // specify which span to process next538 bool F32 = false; // is this an F32 image?539 if (img->type.type == PS_TYPE_F32) {540 F32 = true;541 } else if (img->type.type == PS_TYPE_S32) {542 F32 = false;543 } else { // N.b. You can't trivially add more cases here; F32 is just a bool544 psError(PS_ERR_UNKNOWN, true, "Unsupported psImage type: %d", img->type.type);545 return NULL;546 }547 548 psF32 *imgRowF32 = NULL; // row pointer if F32549 psS32 *imgRowS32 = NULL; // " " " " !F32550 551 const int row0 = img->row0;552 const int col0 = img->col0;553 const int numRows = img->numRows;554 const int numCols = img->numCols;555 556 /********************************************************************************************************/557 558 Old_Spartspan *sspan = NULL;559 for (int i = 0; i < startspans->n; i++) {560 sspan = startspans->data[i];561 if (sspan->direction != PM_SSPAN_DONE) {562 break;563 }564 }565 if (sspan == NULL || sspan->direction == PM_SSPAN_DONE) { // no more Old_Spartspans to process566 return false;567 }568 /*569 * Work570 */571 const PM_SSPAN_DIR dir = sspan->direction;572 /*573 * Set initial span to the startspan574 */575 int x0 = sspan->span->x0 - col0, x1 = sspan->span->x1 - col0;576 /*577 * Go through image identifying objects578 */579 int nx0, nx1 = -1; // new values of x0, x1580 const int di = (dir == PM_SSPAN_UP) ? 1 : -1; // how much i changes to get to the next row581 for (int i = sspan->span->y -row0 + di; i < numRows && i >= 0; i += di) {582 imgRowF32 = img->data.F32[i]; // only one of583 imgRowS32 = img->data.S32[i]; // these is valid!584 //585 // Search left from the pixel diagonally to the left of (i - di, x0). If there's586 // a connected span there it may need to grow up and/or down, so push it onto587 // the stack for later consideration588 //589 nx0 = -1;590 for (int j = x0 - 1; j >= -1; j--) {591 double pixVal = (j < 0) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);592 if (pixVal < threshold) {593 if (j < x0 - 1) { // we found some pixels above threshold594 nx0 = j + 1;595 }596 break;597 }598 }599 600 if (nx0 < 0) { // no span to the left601 nx1 = x0 - 1; // we're going to resume searching at nx1 + 1602 } else {603 //604 // Search right in leftmost span605 //606 //nx1 = 0; // make gcc happy607 for (int j = nx0 + 1; j <= numCols; j++) {608 double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);609 if (pixVal < threshold) {610 nx1 = j - 1;611 break;612 }613 }614 615 const pmSpan *sp = pmFootprintAddSpan(fp, i + row0, nx0 + col0, nx1 + col0);616 617 old_add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART);618 }619 //620 // Now look for spans connected to the old span. The first of these we'll621 // simply process, but others will have to be deferred for later consideration.622 //623 // In fact, if the span overhangs to the right we'll have to defer the overhang624 // until later too, as it too can grow in both directions625 //626 // Note that column numCols exists virtually, and always ends the last span; this627 // is why we claim below that sx1 is always set628 //629 bool first = false; // is this the first new span detected?630 for (int j = nx1 + 1; j <= x1 + 1; j++) {631 double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);632 if (pixVal >= threshold) {633 int sx0 = j++; // span that we're working on is sx0:sx1634 int sx1 = -1; // We know that if we got here, we'll also set sx1635 for (; j <= numCols; j++) {636 double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);637 if (pixVal < threshold) { // end of span638 sx1 = j;639 break;640 }641 }642 assert (sx1 >= 0);643 644 const pmSpan *sp;645 if (first) {646 if (sx1 <= x1) {647 sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);648 old_add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_DONE);649 } else { // overhangs to right650 sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, x1 + col0);651 old_add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_DONE);652 sp = pmFootprintAddSpan(fp, i + row0, x1 + 1 + col0, sx1 + col0 - 1);653 old_add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART);654 }655 first = false;656 } else {657 sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);658 old_add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART);659 }660 }661 }662 663 if (first == false) { // we're done664 break;665 }666 667 x0 = nx0; x1 = nx1;668 }669 /*670 * Cleanup671 */672 673 sspan->direction = PM_SSPAN_DONE;674 return true;675 }676 677 /*678 * Go through an image, starting at (row, col) and assembling all the pixels679 * that are connected to that point (in a chess kings-move sort of way) into680 * a pmFootprint.681 *682 * This is much slower than pmFindFootprints if you want to find lots of683 * footprints, but if you only want a small region about a given point it684 * can be much faster685 *686 * N.b. The returned pmFootprint is not in "normal form"; that is the pmSpans687 * are not sorted by increasing y, x0, x1. If this matters to you, call688 * pmFootprintNormalize()689 */690 pmFootprint *691 pmOldFindFootprintAtPoint(const psImage *img, // image to search692 const float threshold, // Threshold693 int row, int col) { // starting position (in img's parent's coordinate system)694 assert(img != NULL);695 696 bool F32 = false; // is this an F32 image?697 if (img->type.type == PS_TYPE_F32) {698 F32 = true;699 } else if (img->type.type == PS_TYPE_S32) {700 F32 = false;701 } else { // N.b. You can't trivially add more cases here; F32 is just a bool702 psError(PS_ERR_UNKNOWN, true, "Unsupported psImage type: %d", img->type.type);703 return NULL;704 }705 psF32 *imgRowF32 = NULL; // row pointer if F32706 psS32 *imgRowS32 = NULL; // " " " " !F32707 708 const int row0 = img->row0;709 const int col0 = img->col0;710 const int numRows = img->numRows;711 const int numCols = img->numCols;712 /*713 * Is point in image, and above threshold?714 */715 row -= row0; col -= col0;716 if (row < 0 || row >= numRows ||717 col < 0 || col >= numCols) {718 psError(PS_ERR_BAD_PARAMETER_VALUE, true,719 "row/col == (%d, %d) are out of bounds [%d--%d, %d--%d]",720 row + row0, col + col0, row0, row0 + numRows - 1, col0, col0 + numCols - 1);721 return NULL;722 }723 724 double pixVal = F32 ? img->data.F32[row][col] : img->data.S32[row][col];725 if (pixVal < threshold) {726 return pmFootprintAlloc(0, img);727 }728 729 pmFootprint *fp = pmFootprintAlloc(1 + img->numRows/10, img);730 /*731 * Find starting span passing through (row, col)732 */733 psArray *startspans = psArrayAllocEmpty(1); // spans where we have to restart the search734 735 imgRowF32 = img->data.F32[row]; // only one of736 imgRowS32 = img->data.S32[row]; // these is valid!737 {738 int i;739 for (i = col; i >= 0; i--) {740 pixVal = F32 ? imgRowF32[i] : imgRowS32[i];741 if (pixVal < threshold) {742 break;743 }744 }745 int i0 = i;746 for (i = col; i < numCols; i++) {747 pixVal = F32 ? imgRowF32[i] : imgRowS32[i];748 if (pixVal < threshold) {749 break;750 }751 }752 int i1 = i;753 const pmSpan *sp = pmFootprintAddSpan(fp, row + row0, i0 + col0 + 1, i1 + col0 - 1);754 755 old_add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART);756 }757 /*758 * Now workout from those Old_Spartspans, searching for pixels above threshold759 */760 while (old_do_startspan(fp, img, threshold, startspans)) continue;761 /*762 * Cleanup763 */764 psFree(startspans); // restores the image pixel765 766 return fp; // pmFootprint really767 }768 769 /************************************************************************************************************/770 /*771 * A data structure to hold the starting point for a search for pixels above threshold,772 * used by pmFindFootprintAtPoint773 *774 * We don't want to find this span again --- it's already part of the footprint ---775 * so we set the pixels in the image to as small a value as we can. N.b. We may want776 * to reconsider this when dealing with difference images777 */778 #if 0779 typedef enum {PM_SSPAN_DOWN = 0, // scan down from this span780 PM_SSPAN_UP, // scan up from this span781 PM_SSPAN_RESTART, // restart scanning from this span782 PM_SSPAN_DONE // this span is processed783 } PM_SSPAN_DIR; // How to continue searching784 #endif785 786 444 // 787 445 // An enum for mask's pixel values. We're looking for pixels that are above threshold, and 788 446 // we keep extra book-keeping information in the PM_SSPAN_STOP plane. It's simpler to be 789 447 // able to check for 448 // 790 449 enum { 791 450 PM_SSPAN_INITIAL = 0x0, // initial state of pixels. … … 793 452 PM_SSPAN_STOP = 0x2 // you may stop searching when you see this pixel 794 453 }; 795 454 // 455 // The struct that remembers how to [re-]start scanning the image for pixels 456 // 796 457 typedef struct { 797 458 const pmSpan *span; // save the pixel range … … 1440 1101 * starting point, discard the peak. 1441 1102 */ 1442 psErrorCode pmOldFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint1443 const psImage *weight, // corresponding variance image1444 pmFootprint *fp, // Footprint containing mortal peaks1445 const float nsigma_delta, // how many sigma above local background a peak1446 // needs to be to survive1447 const float min_threshold) { // minimum permitted coll height1448 assert (img != NULL); assert (img->type.type == PS_TYPE_F32);1449 assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);1450 assert (img->row0 == weight->row0 && img->col0 == weight->col0);1451 assert (fp != NULL);1452 1453 if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do1454 return PS_ERR_NONE;1455 }1456 1457 psRegion subRegion; // desired subregion; 1 larger than bounding box (grr)1458 subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;1459 subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;1460 const psImage *subImg = psImageSubset((psImage *)img, subRegion);1461 const psImage *subWt = psImageSubset((psImage *)weight, subRegion);1462 assert (subImg != NULL && subWt != NULL);1463 //1464 // The brightest peak is always safe; go through other peaks trying to cull them1465 //1466 for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop1467 const pmPeak *peak = fp->peaks->data[i];1468 int x = peak->x - subImg->col0;1469 int y = peak->y - subImg->row0;1470 //1471 // Find the level nsigma below the peak that must separate the peak1472 // from any of its friends1473 //1474 assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);1475 const float stdev = sqrt(subWt->data.F32[y][x]);1476 float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;1477 if (isnan(threshold) || threshold < min_threshold) {1478 #if 1 // min_threshold is assumed to be below the detection threshold,1479 // so all the peaks are pmFootprint, and this isn't the brightest1480 (void)psArrayRemoveIndex(fp->peaks, i);1481 i--; // we moved everything down one1482 continue;1483 #else1484 #error n.b. We will be running LOTS of checks at this threshold, so only find the footprint once1485 threshold = min_threshold;1486 #endif1487 }1488 if (threshold > subImg->data.F32[y][x]) {1489 threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;1490 }1491 1492 #if 01493 const in npixMin = 1;1494 psArray *peakFootprints = pmFindFootprints(subImg, threshold, npixMin);1495 psImage *idImg = pmSetFootprintArrayIDs(peakFootprints, true);1496 psFree(peakFootprints);1497 const int peak_id = idImg->data.S32[y][x]; // the ID for the peak of interest1498 assert (peak_id != 0);1499 #else1500 const int peak_id = 1; // the ID for the peak of interest1501 pmFootprint *peakFootprint = pmOldFindFootprintAtPoint(subImg, threshold, peak->y, peak->x);1502 psImage *idImg = pmSetFootprintID(peakFootprint, peak_id);1503 psFree(peakFootprint);1504 #endif1505 1506 int j;1507 for (j = 0; j < i; j++) {1508 const pmPeak *peak2 = fp->peaks->data[j];1509 int x2 = peak2->x - subImg->col0;1510 int y2 = peak2->y - subImg->row0;1511 const int peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak1512 1513 if (peak2_id == peak_id) { // There's a brighter peak within the footprint above1514 ; // threshold; so cull our initial peak1515 (void)psArrayRemoveIndex(fp->peaks, i);1516 i--; // we moved everything down one1517 break;1518 }1519 }1520 if (j == i) {1521 j++;1522 }1523 1524 psFree(idImg);1525 }1526 1527 psFree((psImage *)subImg);1528 psFree((psImage *)subWt);1529 1530 return PS_ERR_NONE;1531 }1532 1533 /************************************************************************************************************/1534 /*1535 * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently1536 * isolated. More precisely, for each peak find the highest coll that you'd have to traverse1537 * to reach a still higher peak --- and if that coll's more than nsigma DN below your1538 * starting point, discard the peak.1539 */1540 1103 psErrorCode pmFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint 1541 1104 const psImage *weight, // corresponding variance image … … 1561 1124 // 1562 1125 // We need a psArray of peaks brighter than the current peak. We'll fake this 1563 // for efficiency --- shhh 1126 // by reusing the fp->peaks but lying about n. 1127 // 1128 // We do this for efficiency (otherwise I'd need two peaks lists), and we are 1129 // rather too chummy with psArray in consequence. But it works. 1564 1130 // 1565 1131 psArray *brightPeaks = psArrayAlloc(0); … … 1595 1161 } 1596 1162 1597 const int peak_id = 1; 1163 const int peak_id = 1; // the ID for the peak of interest 1598 1164 brightPeaks->n = i; // only stop at a peak brighter than we are 1599 1165 pmFootprint *peakFootprint = pmFindFootprintAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
Note:
See TracChangeset
for help on using the changeset viewer.
