53#include "boost/format.hpp"
69using IdPixelT = std::uint64_t;
71struct Threshold_traits {};
72struct ThresholdLevel_traits :
public Threshold_traits {
74struct ThresholdPixelLevel_traits :
public Threshold_traits {
76struct ThresholdBitmask_traits :
public Threshold_traits {
79template <
typename PixelT>
82 explicit SetIdImage(std::uint64_t
const id,
bool overwriteId =
false,
long const idMask = 0x0)
85 _withSetReplace(false),
86 _overwriteId(overwriteId),
91 pex::exceptions::InvalidParameterError,
92 str(boost::format(
"Id 0x%x sets bits in the protected mask 0x%x") % _id % _idMask));
96 SetIdImage(std::uint64_t
const id, std::set<std::uint64_t> *oldIds,
bool overwriteId =
false,
97 long const idMask = 0x0)
100 _withSetReplace(true),
101 _overwriteId(overwriteId),
103 _pos(oldIds->
begin()) {
106 pex::exceptions::InvalidParameterError,
107 str(boost::format(
"Id 0x%x sets bits in the protected mask 0x%x") % _id % _idMask));
113 auto val = input & ~_idMask;
115 if (val != 0 && _withSetReplace) {
116 _pos = _oldIds->insert(_pos, val);
119 input = (input & _idMask) + _id;
126 std::uint64_t
const _id;
128 bool _withSetReplace;
130 std::set<std::uint64_t> *_oldIds;
131 std::set<std::uint64_t>::const_iterator _pos;
139inline bool isBadPixel(T) {
144inline bool isBadPixel(
float val) {
149inline bool isBadPixel(
double val) {
156int nbit(
unsigned long i) {
171class FindIdsInFootprint {
173 explicit FindIdsInFootprint() : _ids(), _old(0) {}
189 std::set<T>
const &getIds()
const {
return _ids; }
200 bool operator()(std::shared_ptr<PeakRecord const> a, std::shared_ptr<PeakRecord const> b) {
201 if (a->getPeakValue() !=
b->getPeakValue()) {
202 return (a->getPeakValue() >
b->getPeakValue());
205 if (a->getIx() !=
b->getIx()) {
206 return (a->getIx() <
b->getIx());
209 return (a->getIy() <
b->getIy());
223 bool const circular = ctrl.isCircular().first && ctrl.isCircular().second;
224 bool const isotropic = ctrl.isIsotropic().second;
225 bool const left = ctrl.isLeft().first && ctrl.isLeft().second;
226 bool const right = ctrl.isRight().first && ctrl.isRight().second;
227 bool const up = ctrl.isUp().first && ctrl.isUp().second;
228 bool const down = ctrl.isDown().first && ctrl.isDown().second;
230 lsst::geom::Box2I
const region = lhs.getRegion();
231 if (region != rhs.getRegion()) {
232 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
233 boost::format(
"The two FootprintSets must have the same region").
str());
240 FootprintList
const &lhsFootprints = *lhs.getFootprints();
241 FootprintList
const &rhsFootprints = *rhs.getFootprints();
242 int const nLhs = lhsFootprints.size();
243 int const nRhs = rhsFootprints.size();
245 if (nLhs > 0 && nRhs > 0 &&
246 lhsFootprints[0]->getPeaks().getSchema() != rhsFootprints[0]->getPeaks().getSchema()) {
248 "FootprintSets to be merged must have identical peak schemas.");
255 int const lhsIdNbit = nbit(nLhs);
256 int const lhsIdMask = (lhsIdNbit == 0) ? 0x0 : (1 << lhsIdNbit) - 1;
260 (boost::format(
"%d + %d footprints need too many bits; change IdPixelT typedef") %
269 using OldIdMap = std::map<int, std::set<std::uint64_t>>;
270 OldIdMap overwrittenIds;
272 auto grower = [&circular, &up, &down, &
left, &
right, &isotropic](
273 std::shared_ptr<Footprint>
const &foot,
int amount) -> std::shared_ptr<Footprint> {
275 auto element = isotropic ? geom::Stencil::CIRCLE : geom::Stencil::MANHATTAN;
280 int top = up ? amount : 0;
281 int bottom = down ? amount : 0;
282 int lLimit =
left ? amount : 0;
283 int rLimit =
right ? amount : 0;
285 auto yRange = top + bottom + 1;
286 std::vector<geom::Span> spanList;
289 for (
auto dy = 1; dy <= top; ++dy) {
292 for (
auto dy = -1; dy >= -bottom; --dy) {
296 geom::SpanSet structure(
std::move(spanList));
304 for (FootprintList::const_iterator ptr = lhsFootprints.begin(), end = lhsFootprints.end(); ptr != end;
306 std::shared_ptr<Footprint> foot = *ptr;
308 if (rLhs > 0 && foot->getArea() > 0) {
309 foot = grower(foot, rLhs);
312 std::set<std::uint64_t> overwritten;
314 ->clippedTo(idImage->getBBox())
315 ->applyFunctor(SetIdImage<IdPixelT>(
id, &overwritten,
true), *idImage);
317 if (!overwritten.empty()) {
318 overwrittenIds.insert(overwrittenIds.end(), std::make_pair(id, overwritten));
323 id = (1 << lhsIdNbit);
324 for (FootprintList::const_iterator ptr = rhsFootprints.begin(), end = rhsFootprints.end(); ptr != end;
325 ++ptr,
id += (1 << lhsIdNbit)) {
326 std::shared_ptr<Footprint> foot = *ptr;
328 if (rRhs > 0 && foot->getArea() > 0) {
329 foot = grower(foot, rRhs);
332 std::set<std::uint64_t> overwritten;
334 ->clippedTo(idImage->getBBox())
335 ->applyFunctor(SetIdImage<IdPixelT>(
id, &overwritten,
true, lhsIdMask), *idImage);
337 if (!overwritten.
empty()) {
338 overwrittenIds.insert(overwrittenIds.end(),
std::make_pair(
id, overwritten));
345 schema = lhsFootprints[0]->getPeaks().getSchema();
346 }
else if (nRhs > 0) {
347 schema = rhsFootprints[0]->getPeaks().getSchema();
364 FindIdsInFootprint<IdPixelT> idFinder;
365 for (
const auto& foot : *fs.getFootprints()) {
367 foot->getSpans()->applyFunctor(idFinder, *idImage);
369 std::set<std::uint64_t> lhsFootprintIndxs, rhsFootprintIndxs;
371 for (
unsigned int indx : idFinder.getIds()) {
372 if ((indx & lhsIdMask) > 0) {
373 std::uint64_t i = (indx & lhsIdMask) - 1;
374 lhsFootprintIndxs.
insert(i);
378 OldIdMap::iterator mapPtr = overwrittenIds.find(indx);
379 if (mapPtr != overwrittenIds.end()) {
380 std::set<std::uint64_t> &overwritten = mapPtr->second;
382 for (
unsigned long ptr : overwritten) {
383 lhsFootprintIndxs.
insert((ptr & lhsIdMask) - 1);
390 std::uint64_t i = indx - 1;
391 rhsFootprintIndxs.
insert(i);
395 OldIdMap::iterator mapPtr = overwrittenIds.find(indx);
396 if (mapPtr != overwrittenIds.end()) {
397 std::set<std::uint64_t> &overwritten = mapPtr->second;
399 for (
unsigned long ptr : overwritten) {
400 rhsFootprintIndxs.
insert(ptr - 1);
411 for (
unsigned long i : lhsFootprintIndxs) {
412 assert(i < lhsFootprints.size());
413 PeakCatalog const &oldPeaks = lhsFootprints[i]->getPeaks();
415 int const nold = peaks.size();
416 peaks.insert(peaks.end(), oldPeaks.begin(), oldPeaks.end());
421 peaks.getInternal().end(), SortPeaks());
424 for (
unsigned long i : rhsFootprintIndxs) {
425 assert(i < rhsFootprints.size());
426 PeakCatalog const &oldPeaks = rhsFootprints[i]->getPeaks();
428 int const nold = peaks.size();
429 peaks.insert(peaks.end(), oldPeaks.begin(), oldPeaks.end());
432 peaks.getInternal().end(), SortPeaks());
453struct IdSpanCompare {
454 bool operator()(IdSpan
const &a, IdSpan
const &b)
const {
457 }
else if (a.id >
b.id) {
460 return (a.y <
b.y) ? true :
false;
471 while (
id != aliases[
id]) {
472 resolved =
id = aliases[id];
481template <
typename ImageT>
482void findPeaksInFootprint(ImageT
const &
image,
bool polarity, Footprint &foot,
const long margin = 0) {
483 auto spanSet = foot.getSpans();
484 if (spanSet->size() == 0) {
488 throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
493 const long y_min =
bbox.getMinY() + margin;
494 const long y_max =
bbox.getMaxY() - margin;
495 const long x_min =
bbox.getMinX() + margin;
496 const long x_max =
bbox.getMaxX() - margin;
498 const auto image_x0 =
image.getX0();
499 const auto image_y0 =
image.getY0();
501 for (
auto const &spanIter : *spanSet) {
502 auto y = spanIter.getY() - image_y0;
503 if (((y + image_y0) < y_min) || ((y + image_y0) > y_max)) {
506 for (
auto x = spanIter.getMinX() - image_x0; x <= spanIter.getMaxX() - image_x0; ++x) {
507 if (((x + image_x0) < x_min) || ((x + image_x0) > x_max)) {
510 auto val =
image(x, y);
512 if (
image(x - 1, y + 1) > val ||
image(x, y + 1) > val ||
image(x + 1, y + 1) > val ||
513 image(x - 1, y) > val ||
image(x + 1, y) > val ||
image(x - 1, y - 1) > val ||
514 image(x, y - 1) > val ||
image(x + 1, y - 1) > val) {
518 if (
image(x - 1, y + 1) < val ||
image(x, y + 1) < val ||
image(x + 1, y + 1) < val ||
519 image(x - 1, y) < val ||
image(x + 1, y) < val ||
image(x - 1, y - 1) < val ||
520 image(x, y - 1) < val ||
image(x + 1, y - 1) < val) {
524 foot.addPeak(x + image_x0, y + image_y0, val);
529template <
typename ImageT>
530class FindMaxInFootprint {
532 explicit FindMaxInFootprint(
bool polarity)
533 : _polarity(polarity),
536 _min(std::numeric_limits<double>::
max()),
537 _max(-std::numeric_limits<double>::
max()) {}
555 void addRecord(Footprint &foot)
const { foot.addPeak(_x, _y, _polarity ? _max : _min); }
563template <
typename ImageT,
typename ThresholdT>
564void findPeaks(std::shared_ptr<Footprint> foot, ImageT
const &img,
bool polarity, ThresholdT) {
565 findPeaksInFootprint(img, polarity, *foot, 1);
570 std::stable_sort(foot->getPeaks().getInternal().begin(), foot->getPeaks().getInternal().end(),
573 if (foot->getPeaks().empty()) {
574 FindMaxInFootprint<typename ImageT::Pixel> maxFinder(polarity);
575 foot->getSpans()->applyFunctor(maxFinder,
ndarray::ndImage(img.getArray(), img.getXY0()));
576 maxFinder.addRecord(*foot);
581template <
typename ImageT>
582void findPeaks(std::shared_ptr<Footprint>, ImageT
const &,
bool, ThresholdBitmask_traits) {
590template <
typename ImagePixelT,
typename IterT>
591static inline bool inFootprint(ImagePixelT pixVal, IterT,
bool polarity,
double thresholdVal,
592 ThresholdLevel_traits) {
593 return (polarity ? pixVal : -pixVal) >= thresholdVal;
596template <
typename ImagePixelT,
typename IterT>
597static inline bool inFootprint(ImagePixelT pixVal, IterT var,
bool polarity,
double thresholdVal,
598 ThresholdPixelLevel_traits) {
599 return (polarity ? pixVal : -pixVal) >= thresholdVal * ::sqrt(*var);
602template <
typename ImagePixelT,
typename IterT>
603static inline bool inFootprint(ImagePixelT pixVal, IterT,
bool,
double thresholdVal,
604 ThresholdBitmask_traits) {
605 return (pixVal &
static_cast<long>(thresholdVal));
611template <
typename IterT>
612static inline IterT advancePtr(IterT varPtr, Threshold_traits) {
616template <
typename IterT>
617static inline IterT advancePtr(IterT varPtr, ThresholdPixelLevel_traits) {
625template <
typename ImagePixelT,
typename MaskPixelT,
typename VariancePixelT,
typename ThresholdTraitT>
626static void findFootprints(
627 typename FootprintSet::FootprintList *_footprints,
628 lsst::geom::Box2I
const &_region,
631 double const footprintThreshold,
632 double const includeThresholdMultiplier,
637 PeakTable::makeMinimalSchema()
644 double includeThreshold = footprintThreshold * includeThresholdMultiplier;
646 int const row0 = img.
getY0();
647 int const col0 = img.
getX0();
654 std::vector<int> id1(width + 2);
656 std::vector<int> id2(width + 2);
658 std::vector<int>::iterator idc = id1.begin() + 1;
659 std::vector<int>::iterator idp = id2.begin() + 1;
661 std::vector<int> aliases;
662 aliases.
reserve(1 + height / 20);
664 std::vector<IdSpan> spans;
675 for (
int y = 0;
y != height; ++
y) {
676 if (idc == id1.begin() + 1) {
677 idc = id2.begin() + 1;
678 idp = id1.begin() + 1;
680 idc = id1.begin() + 1;
681 idp = id2.begin() + 1;
686 bool good = (includeThresholdMultiplier == 1.0);
689 x_var_iterator varPtr = (var ==
nullptr) ?
nullptr : var->row_begin(
y);
690 for (
int x = 0;
x < width; ++
x, ++pixPtr, varPtr = advancePtr(varPtr, ThresholdTraitT())) {
691 ImagePixelT
const pixVal = *pixPtr;
693 if (isBadPixel(pixVal) ||
694 !inFootprint(pixVal, varPtr, polarity, footprintThreshold, ThresholdTraitT())) {
702 if (idc[x - 1] != 0) {
704 }
else if (idp[x - 1] != 0) {
706 }
else if (idp[x] != 0) {
708 }
else if (idp[x + 1] != 0) {
723 if (idp[x + 1] != 0 && idp[x + 1] !=
id) {
726 idc[
x] =
id = idp[
x + 1];
729 if (!good && inFootprint(pixVal, varPtr, polarity, includeThreshold, ThresholdTraitT())) {
742 for (
auto & span : spans) {
748 if (spans.size() > 0) {
749 std::sort(spans.begin(), spans.end(), IdSpanCompare());
755 if (spans.size() > 0) {
758 for (
unsigned int i = 0; i <= spans.size(); i++) {
759 if (i == spans.size() || spans[i].id !=
id) {
761 std::vector<geom::Span> tempSpanList;
762 for (; i0 < i; i0++) {
763 good |= spans[i0].good;
764 tempSpanList.
emplace_back(spans[i0].y + row0, spans[i0].x0 + col0, spans[i0].x1 + col0);
769 if (good && fp->getArea() >=
static_cast<std::size_t>(npixMin)) {
770 _footprints->push_back(fp);
774 if (i < spans.size()) {
783 using fiterator = FootprintSet::FootprintList::iterator;
784 for (
auto const &_footprint : *_footprints) {
785 findPeaks(_footprint, img, polarity, ThresholdTraitT());
790template <
typename ImagePixelT>
792 int const npixMin,
bool const setPeaks,
table::Schema const &peakSchema)
793 : _footprints(new FootprintList()), _region(img.getBBox()) {
794 using VariancePixelT =
float;
796 findFootprints<ImagePixelT, image::MaskPixel, VariancePixelT, ThresholdLevel_traits>(
797 _footprints.get(), _region, img,
nullptr, threshold.getValue(img),
798 threshold.getIncludeMultiplier(), threshold.getPolarity(), npixMin, setPeaks, peakSchema);
803template <
typename MaskPixelT>
805 : _footprints(new FootprintList()), _region(msk.getBBox()) {
806 switch (threshold.getType()) {
807 case Threshold::BITMASK:
808 findFootprints<MaskPixelT, MaskPixelT, float, ThresholdBitmask_traits>(
809 _footprints.get(), _region, msk, nullptr, threshold.getValue(),
810 threshold.getIncludeMultiplier(), threshold.getPolarity(), npixMin, false);
813 case Threshold::VALUE:
814 findFootprints<MaskPixelT, MaskPixelT, float, ThresholdLevel_traits>(
815 _footprints.get(), _region, msk, nullptr, threshold.getValue(),
816 threshold.getIncludeMultiplier(), threshold.getPolarity(), npixMin, false);
820 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
821 "You must specify a numerical threshold value with a Mask");
825template <
typename ImagePixelT,
typename MaskPixelT>
827 Threshold
const &threshold,
std::string const &planeName,
int const npixMin,
829 : _footprints(new FootprintList()),
834 switch (threshold.getType()) {
835 case Threshold::PIXEL_STDEV:
836 findFootprints<ImagePixelT, MaskPixelT, VariancePixelT, ThresholdPixelLevel_traits>(
837 _footprints.get(), _region, *maskedImg.getImage(), maskedImg.getVariance().get(),
838 threshold.getValue(maskedImg), threshold.getIncludeMultiplier(), threshold.getPolarity(),
842 findFootprints<ImagePixelT, MaskPixelT, VariancePixelT, ThresholdLevel_traits>(
843 _footprints.get(), _region, *maskedImg.getImage(), maskedImg.getVariance().get(),
844 threshold.getValue(maskedImg), threshold.getIncludeMultiplier(), threshold.getPolarity(),
849 if (planeName ==
"") {
856 mask->addMaskPlane(planeName);
858 MaskPixelT
const bitPlane =
mask->getPlaneBitMask(planeName);
862 for (
auto const &fIter : *_footprints) {
863 fIter->getSpans()->setMask(*(maskedImg.
getMask()), bitPlane);
868 : _footprints(
std::
make_shared<FootprintList>()), _region(region) {}
870FootprintSet::FootprintSet(FootprintSet
const &rhs) : _footprints(new FootprintList), _region(rhs._region) {
871 _footprints->reserve(rhs._footprints->size());
872 for (FootprintSet::FootprintList::const_iterator ptr = rhs._footprints->begin(),
873 end = rhs._footprints->end();
875 _footprints->push_back(std::make_shared<Footprint>(**ptr));
880FootprintSet::FootprintSet(FootprintSet &&rhs) : FootprintSet(rhs) {}
882FootprintSet &FootprintSet::operator=(FootprintSet
const &rhs) {
883 FootprintSet tmp(rhs);
889FootprintSet &FootprintSet::operator=(FootprintSet &&rhs) {
return *
this = rhs; }
891FootprintSet::~FootprintSet() =
default;
893void FootprintSet::merge(FootprintSet
const &rhs,
int tGrow,
int rGrow,
bool isotropic) {
894 FootprintControl
const ctrl(
true, isotropic);
895 FootprintSet fs = mergeFootprintSets(*
this, tGrow, rhs, rGrow, ctrl);
902 for (
auto const &ptr : *_footprints) {
903 ptr->setRegion(region);
907FootprintSet::FootprintSet(FootprintSet
const &rhs,
int r,
bool isotropic)
908 : _footprints(new FootprintList), _region(rhs._region) {
910 FootprintSet fs = rhs;
914 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
915 (boost::format(
"I cannot grow by negative numbers: %d") % r).str());
918 FootprintControl
const ctrl(
true, isotropic);
919 FootprintSet fs = mergeFootprintSets(FootprintSet(rhs.getRegion()), 0, rhs, r, ctrl);
923FootprintSet::FootprintSet(FootprintSet
const &rhs,
int ngrow, FootprintControl
const &ctrl)
924 : _footprints(new FootprintList), _region(rhs._region) {
926 FootprintSet fs = rhs;
929 }
else if (ngrow < 0) {
930 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
931 str(boost::format(
"I cannot grow by negative numbers: %d") % ngrow));
934 FootprintSet fs = mergeFootprintSets(FootprintSet(rhs.getRegion()), 0, rhs, ngrow, ctrl);
938FootprintSet::FootprintSet(FootprintSet
const &fs1, FootprintSet
const &fs2,
bool const)
939 : _footprints(new FootprintList()), _region(fs1._region) {
941 throw LSST_EXCEPT(pex::exceptions::LogicError,
"NOT IMPLEMENTED");
949 for (
auto const &fIter : *_footprints) {
951 fIter->getSpans()->applyFunctor(SetIdImage<FootprintIdPixel>(
id), *im);
957template <
typename ImagePixelT,
typename MaskPixelT>
959 HeavyFootprintCtrl
const *ctrl) {
960 HeavyFootprintCtrl ctrl_s = HeavyFootprintCtrl();
966 for (FootprintList::iterator ptr = _footprints->begin(), end = _footprints->end(); ptr != end; ++ptr) {
967 ptr->reset(
new HeavyFootprint<ImagePixelT, MaskPixelT>(**ptr, mimg, ctrl));
971void FootprintSet::makeSources(afw::table::SourceCatalog &cat)
const {
972 cat.reserve(cat.size() + _footprints->size());
973 for (
auto const &i : *_footprints) {
980 os << rhs.getFootprints()->size() <<
" footprints:\n";
982 for (
const auto &i : *(rhs.getFootprints())) {
983 os << delimiter << *i;
995#define INSTANTIATE(PIXEL) \
996 template FootprintSet::FootprintSet(image::Image<PIXEL> const &, Threshold const &, int const, \
997 bool const, table::Schema const &); \
998 template FootprintSet::FootprintSet(image::MaskedImage<PIXEL, image::MaskPixel> const &, \
999 Threshold const &, std::string const &, int const, bool const); \
1000 template void FootprintSet::makeHeavy(image::MaskedImage<PIXEL, image::MaskPixel> const &, \
1001 HeavyFootprintCtrl const *)
#define INSTANTIATE(FROMSYS, TOSYS)
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
run-length code for part of object
IdSpan(int id, int y, int x0, int x1)
static afw::table::Schema makeMinimalSchema()
Return a minimal schema for Peak tables and records.
A Threshold is used to pass a threshold value to detection algorithms.
The base class for all image classed (Image, Mask, MaskedImage, ...)
int getX0() const
Return the image's column-origin.
int getWidth() const
Return the number of columns in the image.
int getY0() const
Return the image's row-origin.
int getHeight() const
Return the number of rows in the image.
x_iterator row_begin(int y) const
Return an x_iterator to the start of the y'th row.
A class to represent a 2-dimensional array of pixels.
Represent a 2-dimensional array of bitmask pixels.
A class to manipulate images, masks, and variance as a single object.
MaskPtr getMask() const
Return a (shared_ptr to) the MaskedImage's mask.
Defines the fields and offsets for a table.
An integer coordinate rectangle.
int getMinY() const noexcept
void include(Point2I const &point)
Expand this to ensure that this->contains(point).
int getMinX() const noexcept
T emplace_back(T... args)
T inplace_merge(T... args)
std::ostream & operator<<(std::ostream &os, CameraSysPrefix const &detSysPrefix)
int resolve_alias(const std::vector< int > &aliases, int id)
Follow a chain of aliases, returning the final resolved value.
afw::table::CatalogT< PeakRecord > PeakCatalog
std::uint64_t FootprintIdPixel
Pixel type for FootprintSet::insertIntoImage()
Extent< int, 2 > Extent2I
details::ImageNdGetter< T, inA, inB > ndImage(ndarray::Array< T, inA, inB > const &array, lsst::geom::Point2I xy0=lsst::geom::Point2I())
Marks a ndarray to be interpreted as an image when applying a functor from a SpanSet.
decltype(sizeof(void *)) size_t