52 bool comparePrimaryRun(PrimaryRun
const&
first, PrimaryRun
const&
second) {
62 class ComparePrimaryRunY {
64 bool operator()(PrimaryRun
const& pr,
int yval) {
return pr.y < yval; }
65 bool operator()(
int yval, PrimaryRun
const& pr) {
return yval < pr.y; }
68 class ComparePrimaryRunM {
70 bool operator()(PrimaryRun
const& pr,
int mval) {
return pr.m < mval; }
71 bool operator()(
int mval, PrimaryRun
const& pr) {
return mval < pr.m; }
80 bool spansOverlap(Span
const&
a, Span
const&
b,
bool compareY =
true) {
83 yTruth =
a.getY() ==
b.getY();
85 return (yTruth && ((
a.getMaxX() >=
b.getMinX() &&
a.getMinX() <=
b.getMinX()) ||
86 (
b.getMaxX() >=
a.getMinX() &&
b.getMinX() <=
a.getMinX())))
98 bool spansContiguous(Span
const&
a, Span
const&
b,
bool compareY =
true) {
101 yTruth =
a.getY() ==
b.getY();
103 return (yTruth && ((
a.getMaxX() + 1 >=
b.getMinX() &&
a.getMinX() <=
b.getMinX()) ||
104 (
b.getMaxX() + 1 >=
a.getMinX() &&
b.getMinX() <=
a.getMinX())))
120 template <
typename T,
bool invert>
124 auto maskBBox =
mask.getBBox();
125 for (
auto const& spn : spanSet) {
128 bool started =
false;
131 if (
y < maskBBox.getMinY() ||
y > maskBBox.getMaxY()) {
138 int startX =
std::max(spn.getMinX(), maskBBox.getMinX());
139 int endX =
std::min(spn.getMaxX(), maskBBox.getMaxX());
140 for (
int x = startX;
x <= endX; ++
x) {
146 pixelCompare = !pixelCompare;
163 }
else if (started) {
171 return std::make_shared<SpanSet>(
std::move(newVec));
178 SpanSet::SpanSet() : _spanVector(), _bbox(), _area(0) {}
187 for (
int i = beginY; i < _bbox.
getEndY(); ++i) {
188 _spanVector.push_back(
Span(i, beginX, maxX));
195 if (_spanVector.empty()) {
209 if (_spanVector.size() == 0) {
220 void SpanSet::_runNormalize() {
226 std::sort(_spanVector.begin(), _spanVector.end());
231 newSpans.
reserve(_spanVector.size());
233 newSpans.
push_back(*_spanVector.begin());
239 for (
auto iter = ++(_spanVector.begin());
iter != _spanVector.end(); ++
iter) {
240 auto& newSpansEnd = newSpans.
back();
241 if (spansContiguous(newSpansEnd, *
iter)) {
242 newSpansEnd =
Span(newSpansEnd.getY(),
std::min(newSpansEnd.getMinX(),
iter->getMinX()),
253 void SpanSet::_initialize() {
265 int minX = _spanVector[0].getMinX();
266 int maxX = _spanVector[0].getMaxX();
269 for (
const auto& span : _spanVector) {
270 if (span.getMinX() < minX) {
271 minX = span.getMinX();
273 if (span.getMaxX() > maxX) {
274 maxX = span.getMaxX();
277 _area += span.getMaxX() - span.getMinX() + 1;
322 void SpanSet::_label(
325 auto currentIndex = spn.
getY();
326 if (currentIndex > 0) {
328 for (
auto const& tup : sortMap[currentIndex - 1]) {
329 if (!labelVector[tup.first] && spansOverlap(spn, *(tup.second),
false)) {
330 labelVector[tup.first] = currentLabel;
331 _label(*(tup.second), labelVector, currentLabel, sortMap);
335 if (currentIndex <= _spanVector.back().getY() - 1) {
337 for (
auto& tup : sortMap[currentIndex + 1]) {
338 if (!labelVector[tup.first] && spansOverlap(spn, *(tup.second),
false)) {
339 labelVector[tup.first] = currentLabel;
340 _label(*(tup.second), labelVector, currentLabel, sortMap);
353 for (
auto const& spn : _spanVector) {
357 for (
auto const& currentSpan : _spanVector) {
358 if (!labelVector[index]) {
359 labelVector[index] = currentLabel;
360 _label(currentSpan, labelVector, currentLabel, sortMap);
375 auto labeledPair = _makeLabels();
380 return labeledPair.second <= 2;
384 auto labeledPair = _makeLabels();
385 auto labels = labeledPair.first;
386 auto numberOfLabels = labeledPair.second;
394 if (numberOfLabels == 1) {
395 subRegions.
push_back(std::make_shared<SpanSet>());
398 subRegions.
reserve(numberOfLabels - 1);
402 for (
std::size_t i = 0; i < _spanVector.size(); ++i) {
403 subSpanLists[labels[i] - 1].
push_back(_spanVector[i]);
406 for (
std::size_t i = 0; i < numberOfLabels - 1; ++i) {
407 subRegions.
push_back(std::make_shared<SpanSet>(subSpanLists[i]));
422 return makeShift(
x,
y);
427 return makeShift(offset.getX(), offset.getY());
433 tempVec.
reserve(_spanVector.size());
434 for (
auto const& spn : _spanVector) {
437 return std::make_shared<SpanSet>(
std::move(tempVec),
false);
445 for (
auto const& spn : _spanVector) {
452 return std::make_shared<SpanSet>(
std::move(tempVec),
false);
457 for (
auto const& otherSpan :
other) {
458 for (
auto const& spn : _spanVector) {
459 if (spansOverlap(otherSpan, spn)) {
473 for (
auto const& otherSpn :
other) {
475 for (
auto const& spn : _spanVector) {
495 for (
auto& spn : _spanVector) {
506 double xc = 0, yc = 0;
507 for (
auto const& spn : _spanVector) {
508 int const y = spn.
getY();
511 int const npix = x1 - x0 + 1;
514 xc += npix * 0.5 * (x1 + x0);
525 double const xc = cen.getX();
526 double const yc = cen.getY();
528 double sumxx = 0, sumxy = 0, sumyy = 0;
529 for (
auto const& spn : _spanVector) {
530 int const y = spn.
getY();
531 int const x0 = spn.
getX0();
532 int const x1 = spn.
getX1();
533 int const npix = x1 - x0 + 1;
535 for (
int x = x0;
x <= x1; ++
x) {
536 sumxx += (
x - xc) * (
x - xc);
538 sumxy += npix * (0.5 * (x1 + x0) - xc) * (
y - yc);
539 sumyy += npix * (
y - yc) * (
y - yc);
549 return dilated(*stencilToSpanSet);
554 if (
other.size() == 0) {
555 return std::make_shared<SpanSet>(_spanVector.begin(), _spanVector.end(),
false);
561 for (
auto const& spn : _spanVector) {
562 for (
auto const& otherSpn :
other) {
563 int const xmin = spn.
getMinX() + otherSpn.getMinX();
564 int const xmax = spn.
getMaxX() + otherSpn.getMaxX();
565 int const yval = spn.
getY() + otherSpn.getY();
570 return std::make_shared<SpanSet>(
std::move(tempVec));
577 return eroded(*stencilToSpanSet);
582 if (
other.size() == 0 || this->size() == 0) {
583 return std::make_shared<SpanSet>(_spanVector.begin(), _spanVector.end(),
false);
591 for (
auto const& spn : _spanVector) {
593 for (
auto const& otherSpn :
other) {
594 if ((otherSpn.getMaxX() - otherSpn.getMinX()) <= (spn.
getMaxX() - spn.
getMinX())) {
597 int y = spn.
getY() - otherSpn.getY();
608 for (
int y = primaryRuns.
front().y;
y <= primaryRuns.
back().y; ++
y) {
616 auto otherYRange =
other.back().getY() -
other.front().getY() + 1;
617 if (
std::distance(yRange.first, yRange.second) < otherYRange) {
628 for (
int m = 0;
m < otherYRange; ++
m) {
629 auto mRange =
std::equal_range(yRange.first, yRange.second,
m, ComparePrimaryRunM());
630 if ((mRange.first == mRange.second)) {
637 int startX = mRange.first->xmin;
638 int endX = mRange.first->xmax;
639 for (
auto run = mRange.first + 1;
run != mRange.second; ++
run) {
640 if (
run->xmin > endX) {
642 candidateRuns.
push_back(PrimaryRun{
m,
y, startX, endX});
650 candidateRuns.
push_back(PrimaryRun{
m,
y, startX, endX});
659 for (
auto& good : goodRuns) {
660 for (
auto& cand : candidateRuns) {
661 int start =
std::max(good.xmin, cand.xmin);
672 for (
auto&
run : goodRuns) {
676 return std::make_shared<SpanSet>(
std::move(tempVec));
681 return _spanVector ==
other._spanVector;
686 return _spanVector !=
other._spanVector;
695 for (
auto dy = -r; dy <= r; ++dy) {
696 int dx =
static_cast<int>(sqrt(r * r - dy * dy));
697 tempVec.
push_back(
Span(dy + offset.getY(), -dx + offset.getX(), dx + offset.getX()));
701 for (
auto dy = -r; dy <= r; ++dy) {
702 int dx = r -
abs(dy);
703 tempVec.
push_back(
Span(dy + offset.getY(), -dx + offset.getX(), dx + offset.getX()));
707 for (
auto dy = -r; dy <= r; ++dy) {
709 tempVec.
push_back(
Span(dy + offset.getY(), -dx + offset.getX(), dx + offset.getX()));
713 return std::make_shared<SpanSet>(
std::move(tempVec),
false);
718 return std::make_shared<SpanSet>(pr.
begin(), pr.
end());
724 return std::make_shared<SpanSet>();
727 if (
other == *
this) {
728 return std::make_shared<SpanSet>(this->_spanVector);
731 auto otherIter =
other.begin();
732 for (
auto const& spn : _spanVector) {
733 while (otherIter !=
other.end() && otherIter->getY() <= spn.
getY()) {
734 if (spansOverlap(spn, *otherIter)) {
737 auto newSpan =
Span(spn.
getY(), newMin, newMax);
743 return std::make_shared<SpanSet>(
std::move(tempVec));
749 return std::make_shared<SpanSet>(this->
begin(), this->
end());
752 if (other == *
this) {
753 return std::make_shared<SpanSet>();
762 auto otherIter =
other.begin();
763 for (
auto const& spn : _spanVector) {
765 bool spanStarted =
false;
767 while (otherIter !=
other.end() && otherIter->getY() <= spn.
getY()) {
768 if (spansOverlap(spn, *otherIter)) {
784 if (spn.
getMinX() < otherIter->getMinX()) {
791 if (spn.
getMaxX() > otherIter->getMaxX()) {
793 spanBottom = otherIter->getMaxX() + 1;
802 if (spn.
getMaxX() > otherIter->getMaxX()) {
804 spanBottom = otherIter->getMaxX() + 1;
808 if (otherIter->getMaxX() > spn.
getMaxX()) {
825 return std::make_shared<SpanSet>(
std::move(tempVec));
836 tempVec.
insert(tempVec.
end(), _spanVector.begin(), _spanVector.end());
839 return std::make_shared<SpanSet>(
std::move(tempVec));
862 for (
auto const& tc : toPoints) {
877 newBoxPoints.
clear();
882 auto oldBoxPointIter = oldBoxPoints.cbegin();
884 auto p = *oldBoxPointIter;
885 int const xSource =
std::floor(0.5 + p.getX());
886 int const ySource =
std::floor(0.5 + p.getY());
902 return std::make_shared<SpanSet>(
std::move(tempVec));
905 template <
typename ImageT>
914 auto setterFunc = [](
lsst::geom::Point2I const& point, ImageT& out, ImageT in) { out = in; };
918 tmpSpan->applyFunctor(setterFunc,
image,
val);
924 "Footprint Bounds Outside image, set doClip to true");
928 template <
typename T>
931 auto targetArray =
target.getArray();
932 auto xy0 =
target.getBBox().getMin();
935 T bitmask) { maskVal |= bitmask; };
939 template <
typename T>
942 auto targetArray =
target.getArray();
943 auto xy0 =
target.getBBox().getMin();
946 T bitmask) { maskVal &= ~bitmask; };
950 template <
typename T>
952 return maskIntersect<T, false>(*
this,
other, bitmask);
955 template <
typename T>
957 return maskIntersect<T, true>(*
this,
other, bitmask);
960 template <
typename T>
962 auto comparator = [bitmask](T pixelValue) {
return (pixelValue & bitmask); };
964 return union_(*spanSetFromMask);
969 class SpanSetPersistenceHelper {
976 static SpanSetPersistenceHelper
const& get() {
977 static SpanSetPersistenceHelper instance;
982 SpanSetPersistenceHelper(
const SpanSetPersistenceHelper&) =
delete;
983 SpanSetPersistenceHelper& operator=(
const SpanSetPersistenceHelper&) =
delete;
986 SpanSetPersistenceHelper(SpanSetPersistenceHelper&&) =
delete;
987 SpanSetPersistenceHelper& operator=(SpanSetPersistenceHelper&&) =
delete;
990 SpanSetPersistenceHelper()
997 std::string getSpanSetPersistenceName() {
return "SpanSet"; }
999 class SpanSetFactory :
public table::io::PersistableFactory {
1002 CatalogVector
const& catalogs)
const override {
1006 auto spansCatalog = catalogs.front();
1008 auto const&
keys = SpanSetPersistenceHelper::get();
1011 tempVec.
reserve(spansCatalog.size());
1012 for (
auto const&
val : spansCatalog) {
1015 return std::make_shared<SpanSet>(
std::move(tempVec));
1022 SpanSetFactory registration(getSpanSetPersistenceName());
1026 std::string SpanSet::getPersistenceName()
const {
return getSpanSetPersistenceName(); }
1028 void SpanSet::write(OutputArchiveHandle& handle)
const {
1029 auto const&
keys = SpanSetPersistenceHelper::get();
1030 auto spanCat = handle.makeCatalog(
keys.spanSetSchema);
1031 spanCat.reserve(
size());
1032 for (
auto const&
val : *
this) {
1033 auto record = spanCat.addNew();
1034 record->set(
keys.spanY,
val.getY());
1035 record->set(
keys.spanX0,
val.getX0());
1036 record->set(
keys.spanX1,
val.getX1());
1038 handle.saveCatalog(spanCat);
1043 #define INSTANTIATE_IMAGE_TYPE(T) \
1044 template void SpanSet::setImage<T>(image::Image<T> & image, T val, \
1045 lsst::geom::Box2I const& region = lsst::geom::Box2I(), \
1046 bool doClip = false) const;
1048 #define INSTANTIATE_MASK_TYPE(T) \
1049 template void SpanSet::setMask<T>(image::Mask<T> & target, T bitmask) const; \
1050 template void SpanSet::clearMask<T>(image::Mask<T> & target, T bitmask) const; \
1051 template std::shared_ptr<SpanSet> SpanSet::intersect<T>(image::Mask<T> const& other, T bitmask) const; \
1052 template std::shared_ptr<SpanSet> SpanSet::intersectNot<T>(image::Mask<T> const& other, T bitmask) \
1054 template std::shared_ptr<SpanSet> SpanSet::union_<T>(image::Mask<T> const& other, T bitmask) const;