LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
FootprintSet.cc
Go to the documentation of this file.
1// -*- lsst-c++ -*-
2
3/*
4 * LSST Data Management System
5 * Copyright 2008, 2009, 2010 LSST Corporation.
6 *
7 * This product includes software developed by the
8 * LSST Project (http://www.lsst.org/).
9 *
10 * This program is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 3 of the License, or
13 * (at your option) any later version.
14 *
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the LSST License Statement and
21 * the GNU General Public License along with this program. If not,
22 * see <http://www.lsstcorp.org/LegalNotices/>.
23 */
24
25/*
26 * Utilities to detect sets of Footprint%s
27 *
28 * Create and use an lsst::afw::detection::FootprintSet, a collection of pixels above (or below) a threshold
29 * in an Image
30 *
31 * The "collections of pixels" are represented as lsst::afw::detection::Footprint%s, so an example application
32 * would be:
33 *
34 namespace image = lsst::afw::image; namespace detection = lsst::afw::detection;
35
36 image::MaskedImage<float> img(10,20);
37 *img.getImage() = 100;
38
39 detection::FootprintSet<float> sources(img, 10);
40 cout << "Found " << sources.getFootprints()->size() << " sources" << std::endl;
41 */
42#include <cstdint>
43#include <memory>
44#include <algorithm>
45#include <cassert>
46#include <iostream>
47#include <iterator> // for ostream_iterator
48#include <set>
49#include <string>
50#include <map>
51#include <vector>
52
53#include "boost/format.hpp"
54#include "lsst/pex/exceptions.h"
61
62namespace lsst {
63namespace afw {
64namespace detection {
65
66namespace {
68
69using IdPixelT = std::uint64_t; // Type of temporary Images used in merging Footprints
70
71struct Threshold_traits {};
72struct ThresholdLevel_traits : public Threshold_traits { // Threshold is a single number
73};
74struct ThresholdPixelLevel_traits : public Threshold_traits { // Threshold varies from pixel to pixel
75};
76struct ThresholdBitmask_traits : public Threshold_traits { // Threshold ORs with a bitmask
77};
78
79template <typename PixelT>
80class setIdImage {
81public:
82 explicit setIdImage(std::uint64_t const id, bool overwriteId = false, long const idMask = 0x0)
83 : _id(id),
84 _idMask(idMask),
85 _withSetReplace(false),
86 _overwriteId(overwriteId),
87 _oldIds(nullptr),
88 _pos() {
89 if (_id & _idMask) {
90 throw LSST_EXCEPT(
91 pex::exceptions::InvalidParameterError,
92 str(boost::format("Id 0x%x sets bits in the protected mask 0x%x") % _id % _idMask));
93 }
94 }
95
96 setIdImage(std::uint64_t const id, std::set<std::uint64_t> *oldIds, bool overwriteId = false,
97 long const idMask = 0x0)
98 : _id(id),
99 _idMask(idMask),
100 _withSetReplace(true),
101 _overwriteId(overwriteId),
102 _oldIds(oldIds),
103 _pos(oldIds->begin()) {
104 if (_id & _idMask) {
105 throw LSST_EXCEPT(
106 pex::exceptions::InvalidParameterError,
107 str(boost::format("Id 0x%x sets bits in the protected mask 0x%x") % _id % _idMask));
108 }
109 }
110
111 void operator()(lsst::geom::Point2I const &point, PixelT &input) {
112 if (_overwriteId) {
113 auto val = input & ~_idMask;
114
115 if (val != 0 && _withSetReplace) {
116 _pos = _oldIds->insert(_pos, val);
117 }
118
119 input = (input & _idMask) + _id;
120 } else {
121 input += _id;
122 }
123 }
124
125private:
126 std::uint64_t const _id;
127 long const _idMask;
128 bool _withSetReplace;
129 bool _overwriteId;
131 std::set<std::uint64_t>::const_iterator _pos;
132};
133
134//
135// Define our own functions to handle NaN tests; this gives us the
136// option to define a value for e.g. image::MaskPixel or int
137//
138template <typename T>
139inline bool isBadPixel(T) {
140 return false;
141}
142
143template <>
144inline bool isBadPixel(float val) {
145 return std::isnan(val);
146}
147
148template <>
149inline bool isBadPixel(double val) {
150 return std::isnan(val);
151}
152
153/*
154 * Return the number of bits required to represent a unsigned long
155 */
156int nbit(unsigned long i) {
157 int n = 0;
158 while (i > 0) {
159 ++n;
160 i >>= 1;
161 }
162
163 return n;
164}
165/*
166 * Find the list of pixel values that lie in a Footprint
167 *
168 * Used when the Footprints are constructed from an Image containing Footprint indices
169 */
170template <typename T>
171class FindIdsInFootprint {
172public:
173 explicit FindIdsInFootprint() : _ids(), _old(0) {}
174
175 // Reset everything for a new Footprint
176 void reset() {
177 _ids.clear();
178 _old = 0;
179 }
180
181 // Take by copy and not be reference on purpose
182 void operator()(lsst::geom::Point2I const &point, T val) {
183 if (val != _old) {
184 _ids.insert(val);
185 _old = val;
186 }
187 }
188
189 std::set<T> const &getIds() const { return _ids; }
190
191private:
192 std::set<T> _ids;
193 T _old;
194};
195
196/*
197 * Sort peaks by decreasing pixel value. N.b. -ve peaks are sorted the same way as +ve ones
198 */
199struct SortPeaks {
201 if (a->getPeakValue() != b->getPeakValue()) {
202 return (a->getPeakValue() > b->getPeakValue());
203 }
204
205 if (a->getIx() != b->getIx()) {
206 return (a->getIx() < b->getIx());
207 }
208
209 return (a->getIy() < b->getIy());
210 }
211};
212/*
213 * Worker routine for merging two FootprintSets, possibly growing them as we proceed
214 */
215FootprintSet mergeFootprintSets(FootprintSet const &lhs, // the FootprintSet to be merged to
216 int rLhs, // Grow lhs Footprints by this many pixels
217 FootprintSet const &rhs, // the FootprintSet to be merged into lhs
218 int rRhs, // Grow rhs Footprints by this many pixels
219 FootprintControl const &ctrl // Control how the grow is done
220) {
221 using FootprintList = FootprintSet::FootprintList;
222 // The isXXX routines return <isset, value>
223 bool const circular = ctrl.isCircular().first && ctrl.isCircular().second;
224 bool const isotropic = ctrl.isIsotropic().second; // isotropic grow as opposed to a Manhattan metric
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;
229
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());
234 }
235
236 auto idImage = std::make_shared<image::Image<IdPixelT>>(region);
237 idImage->setXY0(region.getMinX(), region.getMinY());
238 *idImage = 0;
239
240 FootprintList const &lhsFootprints = *lhs.getFootprints();
241 FootprintList const &rhsFootprints = *rhs.getFootprints();
242 int const nLhs = lhsFootprints.size();
243 int const nRhs = rhsFootprints.size();
244 // Require equal peak schemas.
245 if (nLhs > 0 && nRhs > 0 &&
246 lhsFootprints[0]->getPeaks().getSchema() != rhsFootprints[0]->getPeaks().getSchema()) {
247 throw LSST_EXCEPT(pex::exceptions::RuntimeError,
248 "FootprintSets to be merged must have identical peak schemas.");
249 }
250
251 /*
252 * In general the lists of Footprints overlap, so we need to make sure that the IDs can be
253 * uniquely recovered from the idImage. We do this by allocating a range of bits to the lhs IDs
254 */
255 int const lhsIdNbit = nbit(nLhs);
256 int const lhsIdMask = (lhsIdNbit == 0) ? 0x0 : (1 << lhsIdNbit) - 1;
257
258 if (std::size_t(nRhs << lhsIdNbit) > std::numeric_limits<IdPixelT>::max() - 1) {
259 throw LSST_EXCEPT(pex::exceptions::OverflowError,
260 (boost::format("%d + %d footprints need too many bits; change IdPixelT typedef") %
261 nLhs % nRhs)
262 .str());
263 }
264 /*
265 * When we insert grown Footprints into the idImage we can potentially overwrite an entire Footprint,
266 * losing any peaks that it might contain. We'll preserve the overwritten Ids in case we need to
267 * get them back (n.b. Footprints that overlap, but both if which survive, will appear in this list)
268 */
270 OldIdMap overwrittenIds; // here's a map from id -> overwritten IDs
271
272 auto grower = [&circular, &up, &down, &left, &right, &isotropic](
273 std::shared_ptr<Footprint> const &foot, int amount) -> std::shared_ptr<Footprint> {
274 if (circular) {
276 auto tmpFoot = std::make_shared<Footprint>(foot->getSpans()->dilated(amount, element),
277 foot->getRegion());
278 return tmpFoot;
279 } else {
280 int top = up ? amount : 0;
281 int bottom = down ? amount : 0;
282 int lLimit = left ? amount : 0;
283 int rLimit = right ? amount : 0;
284
285 auto yRange = top + bottom + 1;
287 spanList.reserve(yRange);
288
289 for (auto dy = 1; dy <= top; ++dy) {
290 spanList.emplace_back(dy, 0, 0);
291 }
292 for (auto dy = -1; dy >= -bottom; --dy) {
293 spanList.emplace_back(dy, 0, 0);
294 }
295 spanList.emplace_back(0, -lLimit, rLimit);
296 geom::SpanSet structure(std::move(spanList));
297 auto tmpFoot =
298 std::make_shared<Footprint>(foot->getSpans()->dilated(structure), foot->getRegion());
299 return tmpFoot;
300 }
301 };
302
303 IdPixelT id = 1; // the ID inserted into the image
304 for (FootprintList::const_iterator ptr = lhsFootprints.begin(), end = lhsFootprints.end(); ptr != end;
305 ++ptr, ++id) {
307
308 if (rLhs > 0 && foot->getArea() > 0) {
309 foot = grower(foot, rLhs);
310 }
311
312 std::set<std::uint64_t> overwritten;
313 foot->getSpans()
314 ->clippedTo(idImage->getBBox())
315 ->applyFunctor(setIdImage<IdPixelT>(id, &overwritten, true), *idImage);
316
317 if (!overwritten.empty()) {
318 overwrittenIds.insert(overwrittenIds.end(), std::make_pair(id, overwritten));
319 }
320 }
321
322 assert(id <= std::size_t(1 << lhsIdNbit));
323 id = (1 << lhsIdNbit);
324 for (FootprintList::const_iterator ptr = rhsFootprints.begin(), end = rhsFootprints.end(); ptr != end;
325 ++ptr, id += (1 << lhsIdNbit)) {
327
328 if (rRhs > 0 && foot->getArea() > 0) {
329 foot = grower(foot, rRhs);
330 }
331
332 std::set<std::uint64_t> overwritten;
333 foot->getSpans()
334 ->clippedTo(idImage->getBBox())
335 ->applyFunctor(setIdImage<IdPixelT>(id, &overwritten, true, lhsIdMask), *idImage);
336
337 if (!overwritten.empty()) {
338 overwrittenIds.insert(overwrittenIds.end(), std::make_pair(id, overwritten));
339 }
340 }
341
342 // Ensure fields added to the input peak schemas are preserved in the output.
343 table::Schema schema;
344 if (nLhs > 0) {
345 schema = lhsFootprints[0]->getPeaks().getSchema();
346 } else if (nRhs > 0) {
347 schema = rhsFootprints[0]->getPeaks().getSchema();
348 } else {
350 }
351 FootprintSet fs(*idImage, Threshold(1), 1, false, schema); // detect all pixels in rhs + lhs
352
353 /*
354 * Now go through the new Footprints looking up and remembering their progenitor's IDs; we'll use
355 * these IDs to merge the peaks in a moment
356 *
357 * We can't do this as we go through the idFinder as the IDs it returns are
358 * (lhsId + 1) | ((rhsId + 1) << nbit)
359 * and, depending on the geometry, values of lhsId and/or rhsId can appear multiple times
360 * (e.g. if nbit is 2, idFinder IDs 0x5 and 0x6 both contain lhsId = 0) so we get duplicates
361 * of peaks. This is not too bad, but it's a bit of a pain to make the lists unique again,
362 * and we avoid this by this two-step process.
363 */
364 FindIdsInFootprint<IdPixelT> idFinder;
365 for (const auto& foot : *fs.getFootprints()) {
366 // find the (mangled) [lr]hsFootprint IDs that contribute to foot
367 foot->getSpans()->applyFunctor(idFinder, *idImage);
368
369 std::set<std::uint64_t> lhsFootprintIndxs, rhsFootprintIndxs; // indexes into [lr]hsFootprints
370
371 for (unsigned int indx : idFinder.getIds()) {
372 if ((indx & lhsIdMask) > 0) {
373 std::uint64_t i = (indx & lhsIdMask) - 1;
374 lhsFootprintIndxs.insert(i);
375 /*
376 * Now allow for Footprints that vanished beneath this one
377 */
378 OldIdMap::iterator mapPtr = overwrittenIds.find(indx);
379 if (mapPtr != overwrittenIds.end()) {
380 std::set<std::uint64_t> &overwritten = mapPtr->second;
381
382 for (unsigned long ptr : overwritten) {
383 lhsFootprintIndxs.insert((ptr & lhsIdMask) - 1);
384 }
385 }
386 }
387 indx >>= lhsIdNbit;
388
389 if (indx > 0) {
390 std::uint64_t i = indx - 1;
391 rhsFootprintIndxs.insert(i);
392 /*
393 * Now allow for Footprints that vanished beneath this one
394 */
395 OldIdMap::iterator mapPtr = overwrittenIds.find(indx);
396 if (mapPtr != overwrittenIds.end()) {
397 std::set<std::uint64_t> &overwritten = mapPtr->second;
398
399 for (unsigned long ptr : overwritten) {
400 rhsFootprintIndxs.insert(ptr - 1);
401 }
402 }
403 }
404 }
405 /*
406 * We now have a complete set of Footprints that contributed to this one, so merge
407 * all their Peaks into the new one
408 */
409 PeakCatalog &peaks = foot->getPeaks();
410
411 for (unsigned long i : lhsFootprintIndxs) {
412 assert(i < lhsFootprints.size());
413 PeakCatalog const &oldPeaks = lhsFootprints[i]->getPeaks();
414
415 int const nold = peaks.size();
416 peaks.insert(peaks.end(), oldPeaks.begin(), oldPeaks.end());
417 // We use getInternal() here to get the vector of shared_ptr that Catalog uses internally,
418 // which causes the STL algorithm to copy pointers instead of PeakRecords (which is what
419 // it'd try to do if we passed Catalog's own iterators).
420 std::inplace_merge(peaks.getInternal().begin(), peaks.getInternal().begin() + nold,
421 peaks.getInternal().end(), SortPeaks());
422 }
423
424 for (unsigned long i : rhsFootprintIndxs) {
425 assert(i < rhsFootprints.size());
426 PeakCatalog const &oldPeaks = rhsFootprints[i]->getPeaks();
427
428 int const nold = peaks.size();
429 peaks.insert(peaks.end(), oldPeaks.begin(), oldPeaks.end());
430 // See note above on why we're using getInternal() here.
431 std::inplace_merge(peaks.getInternal().begin(), peaks.getInternal().begin() + nold,
432 peaks.getInternal().end(), SortPeaks());
433 }
434 idFinder.reset();
435 }
436
437 return fs;
438}
439/*
440 * run-length code for part of object
441 */
442class IdSpan {
443public:
444 explicit IdSpan(int id, int y, int x0, int x1, double good) : id(id), y(y), x0(x0), x1(x1), good(good) {}
445 int id; /* ID for object */
446 int y; /* Row wherein IdSpan dwells */
447 int x0, x1; /* inclusive range of columns */
448 bool good; /* includes a value over the desired threshold? */
449};
450/*
451 * comparison functor; sort by ID then row
452 */
453struct IdSpanCompare {
454 bool operator()(IdSpan const &a, IdSpan const &b) const {
455 if (a.id < b.id) {
456 return true;
457 } else if (a.id > b.id) {
458 return false;
459 } else {
460 return (a.y < b.y) ? true : false;
461 }
462 }
463};
464/*
465 * Follow a chain of aliases, returning the final resolved value.
466 */
467int resolve_alias(std::vector<int> const &aliases, /* list of aliases */
468 int id) { /* alias to look up */
469 int resolved = id; /* resolved alias */
470
471 while (id != aliases[id]) {
472 resolved = id = aliases[id];
473 }
474
475 return (resolved);
476}
478} // namespace
479
480namespace {
481template <typename ImageT>
482void findPeaksInFootprint(ImageT const &image, bool polarity, Footprint &foot, std::size_t const margin = 0) {
483 auto spanSet = foot.getSpans();
484 if (spanSet->size() == 0) {
485 return;
486 }
487 auto bbox = image.getBBox();
488 for (auto const &spanIter : *spanSet) {
489 auto y = spanIter.getY() - image.getY0();
490 if (static_cast<std::size_t>(y + image.getY0()) < bbox.getMinY() + margin ||
491 static_cast<std::size_t>(y + image.getY0()) > bbox.getMaxY() - margin) {
492 continue;
493 }
494 for (auto x = spanIter.getMinX() - image.getX0(); x <= spanIter.getMaxX() - image.getX0(); ++x) {
495 if (static_cast<std::size_t>(x + image.getX0()) < (bbox.getMinX() + margin) ||
496 static_cast<std::size_t>(x + image.getX0()) > (bbox.getMaxX() - margin)) {
497 continue;
498 }
499 auto val = image(x, y);
500 if (polarity) { // look for +ve peaks
501 if (image(x - 1, y + 1) > val || image(x, y + 1) > val || image(x + 1, y + 1) > val ||
502 image(x - 1, y) > val || image(x + 1, y) > val || image(x - 1, y - 1) > val ||
503 image(x, y - 1) > val || image(x + 1, y - 1) > val) {
504 continue;
505 }
506 } else { // look for -ve "peaks" (pits)
507 if (image(x - 1, y + 1) < val || image(x, y + 1) < val || image(x + 1, y + 1) < val ||
508 image(x - 1, y) < val || image(x + 1, y) < val || image(x - 1, y - 1) < val ||
509 image(x, y - 1) < val || image(x + 1, y - 1) < val) {
510 continue;
511 }
512 }
513
514 foot.addPeak(x + image.getX0(), y + image.getY0(), val);
515 }
516 }
517}
518
519template <typename ImageT>
520class FindMaxInFootprint {
521public:
522 explicit FindMaxInFootprint(bool polarity)
523 : _polarity(polarity),
524 _x(0),
525 _y(0),
526 _min(std::numeric_limits<double>::max()),
527 _max(-std::numeric_limits<double>::max()) {}
528
529 void operator()(lsst::geom::Point2I const &point, ImageT const &val) {
530 if (_polarity) {
531 if (val > _max) {
532 _max = val;
533 _x = point.getX();
534 _y = point.getY();
535 }
536 } else {
537 if (val < _min) {
538 _min = val;
539 _x = point.getX();
540 _y = point.getY();
541 }
542 }
543 }
544
545 void addRecord(Footprint &foot) const { foot.addPeak(_x, _y, _polarity ? _max : _min); }
546
547private:
548 bool _polarity;
549 int _x, _y;
550 double _min, _max;
551};
552
553template <typename ImageT, typename ThresholdT>
554void findPeaks(std::shared_ptr<Footprint> foot, ImageT const &img, bool polarity, ThresholdT) {
555 findPeaksInFootprint(img, polarity, *foot, 1);
556
557 // We use getInternal() here to get the vector of shared_ptr that Catalog uses internally,
558 // which causes the STL algorithm to copy pointers instead of PeakRecords (which is what
559 // it'd try to do if we passed Catalog's own iterators).
560 std::stable_sort(foot->getPeaks().getInternal().begin(), foot->getPeaks().getInternal().end(),
561 SortPeaks());
562
563 if (foot->getPeaks().empty()) {
564 FindMaxInFootprint<typename ImageT::Pixel> maxFinder(polarity);
565 foot->getSpans()->applyFunctor(maxFinder, ndarray::ndImage(img.getArray(), img.getXY0()));
566 maxFinder.addRecord(*foot);
567 }
568}
569
570// No need to search for peaks when processing a Mask
571template <typename ImageT>
572void findPeaks(std::shared_ptr<Footprint>, ImageT const &, bool, ThresholdBitmask_traits) {
573 ;
574}
575} // namespace
576
577/*
578 * Functions to determine if a pixel's in a Footprint
579 */
580template <typename ImagePixelT, typename IterT>
581static inline bool inFootprint(ImagePixelT pixVal, IterT, bool polarity, double thresholdVal,
582 ThresholdLevel_traits) {
583 return (polarity ? pixVal : -pixVal) >= thresholdVal;
584}
585
586template <typename ImagePixelT, typename IterT>
587static inline bool inFootprint(ImagePixelT pixVal, IterT var, bool polarity, double thresholdVal,
588 ThresholdPixelLevel_traits) {
589 return (polarity ? pixVal : -pixVal) >= thresholdVal * ::sqrt(*var);
590}
591
592template <typename ImagePixelT, typename IterT>
593static inline bool inFootprint(ImagePixelT pixVal, IterT, bool, double thresholdVal,
594 ThresholdBitmask_traits) {
595 return (pixVal & static_cast<long>(thresholdVal));
596}
597
598/*
599 * Advance the x_iterator to the variance image, when relevant (it may be NULL otherwise)
600 */
601template <typename IterT>
602static inline IterT advancePtr(IterT varPtr, Threshold_traits) {
603 return varPtr;
604}
605
606template <typename IterT>
607static inline IterT advancePtr(IterT varPtr, ThresholdPixelLevel_traits) {
608 return varPtr + 1;
609}
610
611/*
612 * Here's the working routine for the FootprintSet constructors; see documentation
613 * of the constructors themselves
614 */
615template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT, typename ThresholdTraitT>
616static void findFootprints(
617 typename FootprintSet::FootprintList *_footprints, // Footprints
618 lsst::geom::Box2I const &_region, // BBox of pixels that are being searched
619 image::ImageBase<ImagePixelT> const &img, // Image to search for objects
620 image::Image<VariancePixelT> const *var, // img's variance
621 double const footprintThreshold, // threshold value for footprint
622 double const includeThresholdMultiplier, // threshold (relative to footprintThreshold) for inclusion
623 bool const polarity, // if false, search _below_ thresholdVal
624 int const npixMin, // minimum number of pixels in an object
625 bool const setPeaks, // should I set the Peaks list?
626 table::Schema const &peakSchema =
627 PeakTable::makeMinimalSchema() // Schema to use when defining peak catalog.
628) {
629 int id; /* object ID */
630 int in_span; /* object ID of current IdSpan */
631 int nobj = 0; /* number of objects found */
632 int x0 = 0; /* unpacked from a IdSpan */
633
634 double includeThreshold = footprintThreshold * includeThresholdMultiplier; // Threshold for inclusion
635
636 int const row0 = img.getY0();
637 int const col0 = img.getX0();
638 int const height = img.getHeight();
639 int const width = img.getWidth();
640 /*
641 * Storage for arrays that identify objects by ID. We want to be able to
642 * refer to idp[-1] and idp[width], hence the (width + 2)
643 */
644 std::vector<int> id1(width + 2);
645 std::fill(id1.begin(), id1.end(), 0);
646 std::vector<int> id2(width + 2);
647 std::fill(id2.begin(), id2.end(), 0);
648 std::vector<int>::iterator idc = id1.begin() + 1; // object IDs in current/
649 std::vector<int>::iterator idp = id2.begin() + 1; // previous row
650
651 std::vector<int> aliases; // aliases for initially disjoint parts of Footprints
652 aliases.reserve(1 + height / 20); // initial size of aliases
653
654 std::vector<IdSpan> spans; // y:x0,x1 for objects
655 spans.reserve(aliases.capacity()); // initial size of spans
656
657 aliases.push_back(0); // 0 --> 0
658 /*
659 * Go through image identifying objects
660 */
661 using x_iterator = typename image::Image<ImagePixelT>::x_iterator;
662 using x_var_iterator = typename image::Image<VariancePixelT>::x_iterator;
663
664 in_span = 0; // not in a span
665 for (int y = 0; y != height; ++y) {
666 if (idc == id1.begin() + 1) {
667 idc = id2.begin() + 1;
668 idp = id1.begin() + 1;
669 } else {
670 idc = id1.begin() + 1;
671 idp = id2.begin() + 1;
672 }
673 std::fill_n(idc - 1, width + 2, 0);
674
675 in_span = 0; /* not in a span */
676 bool good = (includeThresholdMultiplier == 1.0); /* Span exceeds the threshold? */
677
678 x_iterator pixPtr = img.row_begin(y);
679 x_var_iterator varPtr = (var == nullptr) ? nullptr : var->row_begin(y);
680 for (int x = 0; x < width; ++x, ++pixPtr, varPtr = advancePtr(varPtr, ThresholdTraitT())) {
681 ImagePixelT const pixVal = *pixPtr;
682
683 if (isBadPixel(pixVal) ||
684 !inFootprint(pixVal, varPtr, polarity, footprintThreshold, ThresholdTraitT())) {
685 if (in_span) {
686 spans.emplace_back(in_span, y, x0, x - 1, good);
687
688 in_span = 0;
689 good = false;
690 }
691 } else { /* a pixel to fix */
692 if (idc[x - 1] != 0) {
693 id = idc[x - 1];
694 } else if (idp[x - 1] != 0) {
695 id = idp[x - 1];
696 } else if (idp[x] != 0) {
697 id = idp[x];
698 } else if (idp[x + 1] != 0) {
699 id = idp[x + 1];
700 } else {
701 id = ++nobj;
702 aliases.push_back(id);
703 }
704
705 idc[x] = id;
706 if (!in_span) {
707 x0 = x;
708 in_span = id;
709 }
710 /*
711 * Do we need to merge ID numbers? If so, make suitable entries in aliases[]
712 */
713 if (idp[x + 1] != 0 && idp[x + 1] != id) {
714 aliases[resolve_alias(aliases, idp[x + 1])] = resolve_alias(aliases, id);
715
716 idc[x] = id = idp[x + 1];
717 }
718
719 if (!good && inFootprint(pixVal, varPtr, polarity, includeThreshold, ThresholdTraitT())) {
720 good = true;
721 }
722 }
723 }
724
725 if (in_span) {
726 spans.emplace_back(in_span, y, x0, width - 1, good);
727 }
728 }
729 /*
730 * Resolve aliases; first alias chains, then the IDs in the spans
731 */
732 for (auto & span : spans) {
733 span.id = resolve_alias(aliases, span.id);
734 }
735 /*
736 * Sort spans by ID, so we can sweep through them once
737 */
738 if (spans.size() > 0) {
739 std::sort(spans.begin(), spans.end(), IdSpanCompare());
740 }
741 /*
742 * Build Footprints from spans
743 */
744 unsigned int i0; // initial value of i
745 if (spans.size() > 0) {
746 id = spans[0].id;
747 i0 = 0;
748 for (unsigned int i = 0; i <= spans.size(); i++) { // <= size to catch the last object
749 if (i == spans.size() || spans[i].id != id) {
750 bool good = false; // Span includes pixel sufficient to include footprint in set?
751 std::vector<geom::Span> tempSpanList;
752 for (; i0 < i; i0++) {
753 good |= spans[i0].good;
754 tempSpanList.emplace_back(spans[i0].y + row0, spans[i0].x0 + col0, spans[i0].x1 + col0);
755 }
756 auto tempSpanSet = std::make_shared<geom::SpanSet>(std::move(tempSpanList));
757 auto fp = std::make_shared<Footprint>(tempSpanSet, peakSchema, _region);
758
759 if (good && fp->getArea() >= static_cast<std::size_t>(npixMin)) {
760 _footprints->push_back(fp);
761 }
762 }
763
764 if (i < spans.size()) {
765 id = spans[i].id;
766 }
767 }
768 }
769 /*
770 * Find all peaks within those Footprints
771 */
772 if (setPeaks) {
773 using fiterator = FootprintSet::FootprintList::iterator;
774 for (auto const &_footprint : *_footprints) {
775 findPeaks(_footprint, img, polarity, ThresholdTraitT());
776 }
777 }
778}
779
780template <typename ImagePixelT>
781FootprintSet::FootprintSet(image::Image<ImagePixelT> const &img, Threshold const &threshold,
782 int const npixMin, bool const setPeaks, table::Schema const &peakSchema)
783 : _footprints(new FootprintList()), _region(img.getBBox()) {
784 using VariancePixelT = float;
785
786 findFootprints<ImagePixelT, image::MaskPixel, VariancePixelT, ThresholdLevel_traits>(
787 _footprints.get(), _region, img, nullptr, threshold.getValue(img),
788 threshold.getIncludeMultiplier(), threshold.getPolarity(), npixMin, setPeaks, peakSchema);
789}
790
791// NOTE: not a template to appease swig (see note by instantiations at bottom)
792
793template <typename MaskPixelT>
794FootprintSet::FootprintSet(image::Mask<MaskPixelT> const &msk, Threshold const &threshold, int const npixMin)
795 : _footprints(new FootprintList()), _region(msk.getBBox()) {
796 switch (threshold.getType()) {
797 case Threshold::BITMASK:
798 findFootprints<MaskPixelT, MaskPixelT, float, ThresholdBitmask_traits>(
799 _footprints.get(), _region, msk, nullptr, threshold.getValue(),
800 threshold.getIncludeMultiplier(), threshold.getPolarity(), npixMin, false);
801 break;
802
803 case Threshold::VALUE:
804 findFootprints<MaskPixelT, MaskPixelT, float, ThresholdLevel_traits>(
805 _footprints.get(), _region, msk, nullptr, threshold.getValue(),
806 threshold.getIncludeMultiplier(), threshold.getPolarity(), npixMin, false);
807 break;
808
809 default:
810 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
811 "You must specify a numerical threshold value with a Mask");
812 }
813}
814
815template <typename ImagePixelT, typename MaskPixelT>
816FootprintSet::FootprintSet(const image::MaskedImage<ImagePixelT, MaskPixelT> &maskedImg,
817 Threshold const &threshold, std::string const &planeName, int const npixMin,
818 bool const setPeaks)
819 : _footprints(new FootprintList()),
820 _region(lsst::geom::Point2I(maskedImg.getX0(), maskedImg.getY0()),
821 lsst::geom::Extent2I(maskedImg.getWidth(), maskedImg.getHeight())) {
823 // Find the Footprints
824 switch (threshold.getType()) {
825 case Threshold::PIXEL_STDEV:
826 findFootprints<ImagePixelT, MaskPixelT, VariancePixelT, ThresholdPixelLevel_traits>(
827 _footprints.get(), _region, *maskedImg.getImage(), maskedImg.getVariance().get(),
828 threshold.getValue(maskedImg), threshold.getIncludeMultiplier(), threshold.getPolarity(),
829 npixMin, setPeaks);
830 break;
831 default:
832 findFootprints<ImagePixelT, MaskPixelT, VariancePixelT, ThresholdLevel_traits>(
833 _footprints.get(), _region, *maskedImg.getImage(), maskedImg.getVariance().get(),
834 threshold.getValue(maskedImg), threshold.getIncludeMultiplier(), threshold.getPolarity(),
835 npixMin, setPeaks);
836 break;
837 }
838 // Set Mask if requested
839 if (planeName == "") {
840 return;
841 }
842 //
843 // Define the maskPlane
844 //
846 mask->addMaskPlane(planeName);
847
848 MaskPixelT const bitPlane = mask->getPlaneBitMask(planeName);
849 //
850 // Set the bits where objects are detected
851 //
852 for (auto const &fIter : *_footprints) {
853 fIter->getSpans()->setMask(*(maskedImg.getMask()), bitPlane);
854 }
855}
856
857FootprintSet::FootprintSet(lsst::geom::Box2I region)
858 : _footprints(std::make_shared<FootprintList>()), _region(region) {}
859
860FootprintSet::FootprintSet(FootprintSet const &rhs) : _footprints(new FootprintList), _region(rhs._region) {
861 _footprints->reserve(rhs._footprints->size());
862 for (FootprintSet::FootprintList::const_iterator ptr = rhs._footprints->begin(),
863 end = rhs._footprints->end();
864 ptr != end; ++ptr) {
865 _footprints->push_back(std::make_shared<Footprint>(**ptr));
866 }
867}
868
869// Delegate to copy-constructor for backward-compatibility
870FootprintSet::FootprintSet(FootprintSet &&rhs) : FootprintSet(rhs) {}
871
872FootprintSet &FootprintSet::operator=(FootprintSet const &rhs) {
873 FootprintSet tmp(rhs);
874 swap(tmp); // See Meyers, Effective C++, Item 11
875 return *this;
876}
877
878// Delegate to copy-assignment for backward-compatibility
879FootprintSet &FootprintSet::operator=(FootprintSet &&rhs) { return *this = rhs; }
880
881FootprintSet::~FootprintSet() = default;
882
883void FootprintSet::merge(FootprintSet const &rhs, int tGrow, int rGrow, bool isotropic) {
884 FootprintControl const ctrl(true, isotropic);
885 FootprintSet fs = mergeFootprintSets(*this, tGrow, rhs, rGrow, ctrl);
886 swap(fs); // Swap the new FootprintSet into place
887}
888
889void FootprintSet::setRegion(lsst::geom::Box2I const &region) {
890 _region = region;
891
892 for (auto const &ptr : *_footprints) {
893 ptr->setRegion(region);
894 }
895}
896
897FootprintSet::FootprintSet(FootprintSet const &rhs, int r, bool isotropic)
898 : _footprints(new FootprintList), _region(rhs._region) {
899 if (r == 0) {
900 FootprintSet fs = rhs;
901 swap(fs); // Swap the new FootprintSet into place
902 return;
903 } else if (r < 0) {
904 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
905 (boost::format("I cannot grow by negative numbers: %d") % r).str());
906 }
907
908 FootprintControl const ctrl(true, isotropic);
909 FootprintSet fs = mergeFootprintSets(FootprintSet(rhs.getRegion()), 0, rhs, r, ctrl);
910 swap(fs); // Swap the new FootprintSet into place
911}
912
913FootprintSet::FootprintSet(FootprintSet const &rhs, int ngrow, FootprintControl const &ctrl)
914 : _footprints(new FootprintList), _region(rhs._region) {
915 if (ngrow == 0) {
916 FootprintSet fs = rhs;
917 swap(fs); // Swap the new FootprintSet into place
918 return;
919 } else if (ngrow < 0) {
920 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
921 str(boost::format("I cannot grow by negative numbers: %d") % ngrow));
922 }
923
924 FootprintSet fs = mergeFootprintSets(FootprintSet(rhs.getRegion()), 0, rhs, ngrow, ctrl);
925 swap(fs); // Swap the new FootprintSet into place
926}
927
928FootprintSet::FootprintSet(FootprintSet const &fs1, FootprintSet const &fs2, bool const)
929 : _footprints(new FootprintList()), _region(fs1._region) {
930 _region.include(fs2._region);
931 throw LSST_EXCEPT(pex::exceptions::LogicError, "NOT IMPLEMENTED");
932}
933
934std::shared_ptr<image::Image<FootprintIdPixel>> FootprintSet::insertIntoImage() const {
935 auto im = std::make_shared<image::Image<FootprintIdPixel>>(_region);
936 *im = 0;
937
938 FootprintIdPixel id = 0;
939 for (auto const &fIter : *_footprints) {
940 id++;
941 fIter->getSpans()->applyFunctor(setIdImage<FootprintIdPixel>(id), *im);
942 }
943
944 return im;
945}
946
947template <typename ImagePixelT, typename MaskPixelT>
948void FootprintSet::makeHeavy(image::MaskedImage<ImagePixelT, MaskPixelT> const &mimg,
949 HeavyFootprintCtrl const *ctrl) {
950 HeavyFootprintCtrl ctrl_s = HeavyFootprintCtrl();
951
952 if (!ctrl) {
953 ctrl = &ctrl_s;
954 }
955
956 for (FootprintList::iterator ptr = _footprints->begin(), end = _footprints->end(); ptr != end; ++ptr) {
957 ptr->reset(new HeavyFootprint<ImagePixelT, MaskPixelT>(**ptr, mimg, ctrl));
958 }
959}
960
961void FootprintSet::makeSources(afw::table::SourceCatalog &cat) const {
962 for (auto const &i : *_footprints) {
964 r->setFootprint(i);
965 }
966}
967
968std::ostream &operator<<(std::ostream &os, FootprintSet const &rhs) {
969 os << rhs.getFootprints()->size() << " footprints:\n";
970 auto delimiter = "";
971 for (const auto &i : *(rhs.getFootprints())) {
972 os << delimiter << *i;
973 delimiter = "\n";
974 }
975 return os;
976}
977
978//
979// Explicit instantiations
980//
981
982#ifndef DOXYGEN
983
984#define INSTANTIATE(PIXEL) \
985 template FootprintSet::FootprintSet(image::Image<PIXEL> const &, Threshold const &, int const, \
986 bool const, table::Schema const &); \
987 template FootprintSet::FootprintSet(image::MaskedImage<PIXEL, image::MaskPixel> const &, \
988 Threshold const &, std::string const &, int const, bool const); \
989 template void FootprintSet::makeHeavy(image::MaskedImage<PIXEL, image::MaskPixel> const &, \
990 HeavyFootprintCtrl const *)
991
992template FootprintSet::FootprintSet(image::Mask<image::MaskPixel> const &, Threshold const &, int const);
993
994template void FootprintSet::setMask(image::Mask<image::MaskPixel> *, std::string const &);
995template void FootprintSet::setMask(std::shared_ptr<image::Mask<image::MaskPixel>>, std::string const &);
996
998INSTANTIATE(int);
999INSTANTIATE(float);
1000INSTANTIATE(double);
1001} // namespace detection
1002} // namespace afw
1003} // namespace lsst
1004
1005#endif // !DOXYGEN
AmpInfoBoxKey bbox
Definition Amplifier.cc:117
double element[2]
Definition BaseTable.cc:90
int end
int max
table::Key< int > id
Definition Detector.cc:162
#define INSTANTIATE(FROMSYS, TOSYS)
Definition Detector.cc:509
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
afw::table::Key< afw::table::Array< MaskPixelT > > mask
uint64_t * ptr
Definition RangeSet.cc:95
std::ostream * os
Definition Schema.cc:557
int y
Definition SpanSet.cc:48
table::Key< int > b
table::Schema schema
Definition python.h:134
T begin(T... args)
T capacity(T... args)
FootprintSet(image::Image< ImagePixelT > const &img, Threshold const &threshold, int const npixMin=1, bool const setPeaks=true, table::Schema const &peakSchema=PeakTable::makeMinimalSchema())
Find a FootprintSet given an Image and a threshold.
std::vector< std::shared_ptr< Footprint > > FootprintList
The FootprintSet's set of Footprints.
IdSpan(int id, int y, int x0, int x1)
Definition CR.cc:68
static afw::table::Schema makeMinimalSchema()
Return a minimal schema for Peak tables and records.
Definition Peak.h:137
The base class for all image classed (Image, Mask, MaskedImage, ...)
Definition ImageBase.h:102
int getX0() const
Return the image's column-origin.
Definition ImageBase.h:306
int getWidth() const
Return the number of columns in the image.
Definition ImageBase.h:294
int getY0() const
Return the image's row-origin.
Definition ImageBase.h:314
int getHeight() const
Return the number of rows in the image.
Definition ImageBase.h:296
x_iterator row_begin(int y) const
Return an x_iterator to the start of the y'th row.
Definition ImageBase.h:401
A class to represent a 2-dimensional array of pixels.
Definition Image.h:51
Represent a 2-dimensional array of bitmask pixels.
Definition Mask.h:77
A class to manipulate images, masks, and variance as a single object.
Definition MaskedImage.h:74
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage's variance.
MaskPtr getMask() const
Return a (shared_ptr to) the MaskedImage's mask.
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
An integer coordinate rectangle.
Definition Box.h:55
int getMinY() const noexcept
Definition Box.h:158
void include(Point2I const &point)
Expand this to ensure that this->contains(point).
Definition Box.cc:152
int getMinX() const noexcept
Definition Box.h:157
T emplace_back(T... args)
T empty(T... args)
T end(T... args)
T fill(T... args)
T fill_n(T... args)
T get(T... args)
T inplace_merge(T... args)
T insert(T... args)
T isnan(T... args)
T left(T... args)
T make_pair(T... args)
T make_shared(T... args)
T move(T... args)
int resolve_alias(const std::vector< int > &aliases, int id)
Follow a chain of aliases, returning the final resolved value.
Definition CR.cc:96
std::uint64_t FootprintIdPixel
Pixel type for FootprintSet::insertIntoImage()
Extent< int, 2 > Extent2I
Definition Extent.h:397
Point< int, 2 > Point2I
Definition Point.h:321
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.
STL namespace.
T push_back(T... args)
T reserve(T... args)
T size(T... args)
T sort(T... args)
ImageT val
Definition CR.cc:146
T stable_sort(T... args)