LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
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 
45 #include <algorithm>
46 #include <cassert>
47 #include <set>
48 #include <string>
49 #include <typeinfo>
50 #include "boost/format.hpp"
51 #include "lsst/pex/exceptions.h"
52 #include "lsst/pex/logging/Trace.h"
60 
61 namespace detection = lsst::afw::detection;
62 namespace image = lsst::afw::image;
63 namespace math = lsst::afw::math;
64 namespace pexLogging = lsst::pex::logging;
65 namespace geom = lsst::afw::geom;
66 
67 /************************************************************************************************************/
68 namespace {
70 
71  typedef boost::uint64_t IdPixelT; // Type of temporary Images used in merging Footprints
72 
73  struct Threshold_traits {
74  };
75  struct ThresholdLevel_traits : public Threshold_traits { // Threshold is a single number
76  };
77  struct ThresholdPixelLevel_traits : public Threshold_traits { // Threshold varies from pixel to pixel
78  };
79  struct ThresholdBitmask_traits : public Threshold_traits { // Threshold ORs with a bitmask
80  };
81 
82  //
83  // Define our own functions to handle NaN tests; this gives us the
84  // option to define a value for e.g. image::MaskPixel or int
85  //
86  template<typename T>
87  inline bool isBadPixel(T) {
88  return false;
89  }
90 
91  template<>
92  inline bool isBadPixel(float val) {
93  return std::isnan(val);
94  }
95 
96  template<>
97  inline bool isBadPixel(double val) {
98  return std::isnan(val);
99  }
100 
101  /*
102  * Return the number of bits required to represent a unsigned long
103  */
104  int nbit(unsigned long i) {
105  int n = 0;
106  while (i > 0) {
107  ++n;
108  i >>= 1;
109  }
110 
111  return n;
112  }
113  /*
114  * Find the list of pixel values that lie in a Footprint
115  *
116  * Used when the Footprints are constructed from an Image containing Footprint indices
117  */
118  template <typename ImageT>
119  class FindIdsInFootprint: public detection::FootprintFunctor<ImageT> {
120  public:
121  explicit FindIdsInFootprint(ImageT const& image
122  ) : detection::FootprintFunctor<ImageT>(image), _ids(), _old(0) {}
124  void reset() {
125  _ids.clear();
126  _old = 0;
127  }
128 
130  void operator()(typename ImageT::xy_locator loc,
131  int x,
132  int y
133  ) {
134  typename ImageT::Pixel val = loc(0, 0);
135 
136  if (val != _old) {
137  _ids.insert(val);
138  _old = val;
139  }
140  }
141 
142  std::set<typename ImageT::Pixel> const& getIds() const {
143  return _ids;
144  }
145 
146  private:
147  std::set<typename ImageT::Pixel> _ids;
148  typename ImageT::Pixel _old;
149  };
150 
151  /********************************************************************************************************/
152  /*
153  * Sort peaks by decreasing pixel value. N.b. -ve peaks are sorted the same way as +ve ones
154  */
155  struct SortPeaks {
157  if (a->getPeakValue() != b->getPeakValue()) {
158  return (a->getPeakValue() > b->getPeakValue());
159  }
160 
161  if (a->getIx() != b->getIx()) {
162  return (a->getIx() < b->getIx());
163  }
164 
165  return (a->getIy() < b->getIy());
166  }
167  };
168  /********************************************************************************************************/
169  /*
170  * Worker routine for merging two FootprintSets, possibly growing them as we proceed
171  */
173  mergeFootprintSets(
174  detection::FootprintSet const &lhs, // the FootprintSet to be merged to
175  int rLhs, // Grow lhs Footprints by this many pixels
176  detection::FootprintSet const &rhs, // the FootprintSet to be merged into lhs
177  int rRhs, // Grow rhs Footprints by this many pixels
178  detection::FootprintControl const& ctrl // Control how the grow is done
179  )
180  {
183  // The isXXX routines return <isset, value>
184  bool const circular = ctrl.isCircular().first && ctrl.isCircular().second;
185  bool const isotropic = ctrl.isIsotropic().second; // isotropic grow as opposed to a Manhattan metric
186  // n.b. Isotropic grows are significantly slower
187  bool const left = ctrl.isLeft().first && ctrl.isLeft().second;
188  bool const right = ctrl.isRight().first && ctrl.isRight().second;
189  bool const up = ctrl.isUp().first && ctrl.isUp().second;
190  bool const down = ctrl.isDown().first && ctrl.isDown().second;
191 
192  geom::Box2I const region = lhs.getRegion();
193  if (region != rhs.getRegion()) {
194  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
195  boost::format("The two FootprintSets must have the same region").str());
196  }
197 
199  idImage->setXY0(region.getMinX(), region.getMinY());
200  *idImage = 0;
201 
202  FootprintList const& lhsFootprints = *lhs.getFootprints();
203  FootprintList const& rhsFootprints = *rhs.getFootprints();
204  int const nLhs = lhsFootprints.size();
205  int const nRhs = rhsFootprints.size();
206  /*
207  * In general the lists of Footprints overlap, so we need to make sure that the IDs can be
208  * uniquely recovered from the idImage. We do this by allocating a range of bits to the lhs IDs
209  */
210  int const lhsIdNbit = nbit(nLhs);
211  int const lhsIdMask = (lhsIdNbit == 0) ? 0x0 : (1 << lhsIdNbit) - 1;
212 
213  if (std::size_t(nRhs << lhsIdNbit) > std::numeric_limits<IdPixelT>::max() - 1) {
214  throw LSST_EXCEPT(lsst::pex::exceptions::OverflowError,
215  (boost::format("%d + %d footprints need too many bits; change IdPixelT typedef")
216  % nLhs % nRhs).str());
217  }
218  /*
219  * When we insert grown Footprints into the idImage we can potentially overwrite an entire Footprint,
220  * losing any peaks that it might contain. We'll preserve the overwritten Ids in case we need to
221  * get them back (n.b. Footprints that overlap, but both if which survive, will appear in this list)
222  */
223  typedef std::map<int, std::set<boost::uint64_t> > OldIdMap;
224  OldIdMap overwrittenIds; // here's a map from id -> overwritten IDs
225 
226  IdPixelT id = 1; // the ID inserted into the image
227  for (FootprintList::const_iterator ptr = lhsFootprints.begin(), end = lhsFootprints.end();
228  ptr != end; ++ptr, ++id) {
229  CONST_PTR(Footprint) foot = *ptr;
230 
231  if (rLhs > 0) {
232  foot = circular ?
233  growFootprint(*foot, rLhs, isotropic) : growFootprint(*foot, rLhs, left, right, up, down);
234  }
235 
236  std::set<boost::uint64_t> overwritten;
237  foot->insertIntoImage(*idImage, id, true, 0x0, &overwritten);
238 
239  if (!overwritten.empty()) {
240  overwrittenIds.insert(overwrittenIds.end(), std::make_pair(id, overwritten));
241  }
242  }
243 
244  assert (id <= std::size_t(1 << lhsIdNbit));
245  id = (1 << lhsIdNbit);
246  for (FootprintList::const_iterator ptr = rhsFootprints.begin(), end = rhsFootprints.end();
247  ptr != end; ++ptr, id += (1 << lhsIdNbit)) {
248  CONST_PTR(Footprint) foot = *ptr;
249 
250  if (rRhs > 0) {
251  foot = circular ?
252  growFootprint(*foot, rRhs, isotropic) : growFootprint(*foot, rRhs, left, right, up, down);
253  }
254 
255  std::set<boost::uint64_t> overwritten;
256  foot->insertIntoImage(*idImage, id, true, lhsIdMask, &overwritten);
257 
258  if (!overwritten.empty()) {
259  overwrittenIds.insert(overwrittenIds.end(), std::make_pair(id, overwritten));
260  }
261  }
262 
264  1, false); // detect all pixels in rhs + lhs
265  /*
266  * Now go through the new Footprints looking up and remembering their progenitor's IDs; we'll use
267  * these IDs to merge the peaks in a moment
268  *
269  * We can't do this as we go through the idFinder as the IDs it returns are
270  * (lhsId + 1) | ((rhsId + 1) << nbit)
271  * and, depending on the geometry, values of lhsId and/or rhsId can appear multiple times
272  * (e.g. if nbit is 2, idFinder IDs 0x5 and 0x6 both contain lhsId = 0) so we get duplicates
273  * of peaks. This is not too bad, but it's a bit of a pain to make the lists unique again,
274  * and we avoid this by this two-step process.
275  */
276  FindIdsInFootprint<image::Image<IdPixelT> > idFinder(*idImage);
277  for (FootprintList::iterator ptr = fs.getFootprints()->begin(),
278  end = fs.getFootprints()->end(); ptr != end; ++ptr) {
279  PTR(Footprint) foot = *ptr;
280 
281  idFinder.apply(*foot); // find the (mangled) [lr]hsFootprint IDs that contribute to foot
282 
283  std::set<boost::uint64_t> lhsFootprintIndxs, rhsFootprintIndxs; // indexes into [lr]hsFootprints
284 
285  for (std::set<IdPixelT>::iterator idptr = idFinder.getIds().begin(),
286  idend = idFinder.getIds().end(); idptr != idend; ++idptr) {
287  unsigned int indx = *idptr;
288  if ((indx & lhsIdMask) > 0) {
289  boost::uint64_t i = (indx & lhsIdMask) - 1;
290  lhsFootprintIndxs.insert(i);
291  /*
292  * Now allow for Footprints that vanished beneath this one
293  */
294  OldIdMap::iterator mapPtr = overwrittenIds.find(indx);
295  if (mapPtr != overwrittenIds.end()) {
296  std::set<boost::uint64_t> &overwritten = mapPtr->second;
297 
298  for (std::set<boost::uint64_t>::iterator ptr = overwritten.begin(),
299  end = overwritten.end(); ptr != end; ++ptr){
300  lhsFootprintIndxs.insert((*ptr & lhsIdMask) - 1);
301  }
302  }
303  }
304  indx >>= lhsIdNbit;
305 
306  if (indx > 0) {
307  boost::uint64_t i = indx - 1;
308  rhsFootprintIndxs.insert(i);
309  /*
310  * Now allow for Footprints that vanished beneath this one
311  */
312  OldIdMap::iterator mapPtr = overwrittenIds.find(indx);
313  if (mapPtr != overwrittenIds.end()) {
314  std::set<boost::uint64_t> &overwritten = mapPtr->second;
315 
316  for (std::set<boost::uint64_t>::iterator ptr = overwritten.begin(),
317  end = overwritten.end(); ptr != end; ++ptr) {
318  rhsFootprintIndxs.insert(*ptr - 1);
319  }
320  }
321  }
322  }
323  /*
324  * We now have a complete set of Footprints that contributed to this one, so merge
325  * all their Peaks into the new one
326  */
327  detection::PeakCatalog &peaks = foot->getPeaks();
328 
329  for (std::set<boost::uint64_t>::iterator ptr = lhsFootprintIndxs.begin(),
330  end = lhsFootprintIndxs.end(); ptr != end; ++ptr) {
331  boost::uint64_t i = *ptr;
332  assert (i < lhsFootprints.size());
333  detection::PeakCatalog const& oldPeaks = lhsFootprints[i]->getPeaks();
334 
335  int const nold = peaks.size();
336  peaks.insert(peaks.end(), oldPeaks.begin(), oldPeaks.end());
337  // We use getInternal() here to get the vector of shared_ptr that Catalog uses internally,
338  // which causes the STL algorithm to copy pointers instead of PeakRecords (which is what
339  // it'd try to do if we passed Catalog's own iterators).
340  std::inplace_merge(peaks.getInternal().begin(), peaks.getInternal().begin() + nold,
341  peaks.getInternal().end(), SortPeaks());
342  }
343 
344  for (std::set<boost::uint64_t>::iterator ptr = rhsFootprintIndxs.begin(),
345  end = rhsFootprintIndxs.end(); ptr != end; ++ptr) {
346  boost::uint64_t i = *ptr;
347  assert (i < rhsFootprints.size());
348  detection::PeakCatalog const& oldPeaks = rhsFootprints[i]->getPeaks();
349 
350  int const nold = peaks.size();
351  peaks.insert(peaks.end(), oldPeaks.begin(), oldPeaks.end());
352  // See note above on why we're using getInternal() here.
353  std::inplace_merge(peaks.getInternal().begin(), peaks.getInternal().begin() + nold,
354  peaks.getInternal().end(), SortPeaks());
355  }
356  }
357 
358  return fs;
359  }
360 /*
361  * run-length code for part of object
362  */
363  class IdSpan {
364  public:
365  typedef boost::shared_ptr<IdSpan> Ptr;
366 
367  explicit IdSpan(int id, int y, int x0, int x1, double good) :
368  id(id), y(y), x0(x0), x1(x1), good(good) {}
369  int id; /* ID for object */
370  int y; /* Row wherein IdSpan dwells */
371  int x0, x1; /* inclusive range of columns */
372  bool good; /* includes a value over the desired threshold? */
373  };
374 /*
375  * comparison functor; sort by ID then row
376  */
377  struct IdSpanCompar : public std::binary_function<const IdSpan::Ptr, const IdSpan::Ptr, bool> {
378  bool operator()(IdSpan::Ptr const a, IdSpan::Ptr const b) {
379  if (a->id < b->id) {
380  return true;
381  } else if (a->id > b->id) {
382  return false;
383  } else {
384  return (a->y < b->y) ? true : false;
385  }
386  }
387  };
388 /*
389  * Follow a chain of aliases, returning the final resolved value.
390  */
391  int resolve_alias(std::vector<int> const &aliases, /* list of aliases */
392  int id) { /* alias to look up */
393  int resolved = id; /* resolved alias */
394 
395  while (id != aliases[id]) {
396  resolved = id = aliases[id];
397  }
398 
399  return(resolved);
400  }
402 }
403 
404 /************************************************************************************************************/
405 
406 namespace {
407  /*
408  * Find all the Peaks within a Footprint
409  */
410  template <typename ImageT>
411  class FindPeaksInFootprint: public detection::FootprintFunctor<ImageT> {
412  public:
413  explicit FindPeaksInFootprint(ImageT const& image,
414  bool polarity,
416  ) : detection::FootprintFunctor<ImageT>(image),
417  _polarity(polarity), _peaks(peaks) {}
418 
420  void operator()(typename ImageT::xy_locator loc,
421  int x,
422  int y
423  ) {
424  typename ImageT::Pixel val = loc(0, 0);
425 
426  if (_polarity) { // look for +ve peaks
427  if (loc(-1, 1) > val || loc( 0, 1) > val || loc( 1, 1) > val ||
428  loc(-1, 0) > val || loc( 1, 0) > val ||
429  loc(-1, -1) > val || loc( 0, -1) > val || loc( 1, -1) > val) {
430  return;
431  }
432  } else { // look for -ve "peaks" (pits)
433  if (loc(-1, 1) < val || loc( 0, 1) < val || loc( 1, 1) < val ||
434  loc(-1, 0) < val || loc( 1, 0) < val ||
435  loc(-1, -1) < val || loc( 0, -1) < val || loc( 1, -1) < val) {
436  return;
437  }
438  }
439 
440  PTR(detection::PeakRecord) newPeak = _peaks.addNew();
441  newPeak->setIx(x);
442  newPeak->setIy(y);
443  newPeak->setFx(x);
444  newPeak->setFy(y);
445  newPeak->setPeakValue(val);
446  }
447  private:
448  bool _polarity;
449  detection::PeakCatalog &_peaks;
450  };
451 
452  /*
453  * Find the maximum (or minimum, if polarity is false) pixel in a Footprint
454  */
455  template <typename ImageT>
456  class FindMaxInFootprint : public detection::FootprintFunctor<ImageT> {
457  public:
458  explicit FindMaxInFootprint(ImageT const& image,
459  bool polarity
460  ) : detection::FootprintFunctor<ImageT>(image),
461  _polarity(polarity), _x(0), _y(0),
462  _min( std::numeric_limits<double>::max()),
463  _max(-std::numeric_limits<double>::max()) {}
464 
466  void reset() {
467  _x = _y = 0;
468  _min = std::numeric_limits<double>::max();
469  _max = -std::numeric_limits<double>::max();
470  }
471  virtual void reset(detection::Footprint const&) {}
472 
474  void operator()(typename ImageT::xy_locator loc,
475  int x,
476  int y
477  ) {
478  typename ImageT::Pixel val = loc(0, 0);
479 
480  if (_polarity) {
481  if (val > _max) {
482  _max = val;
483  _x = x;
484  _y = y;
485  }
486  } else {
487  if (val < _min) {
488  _min = val;
489  _x = x;
490  _y = y;
491  }
492  }
493  }
494 
495  // Add the Footprint's peak to the given PeakCatalog
496  void addPeak(detection::PeakCatalog & peakCat) const {
497  PTR(detection::PeakRecord) newPeak = peakCat.addNew();
498  newPeak->setIx(_x);
499  newPeak->setIy(_y);
500  newPeak->setFx(_x);
501  newPeak->setFy(_y);
502  newPeak->setPeakValue(_polarity ? _max : _min);
503  }
504  private:
505  bool _polarity;
506  int _x, _y;
507  double _min, _max;
508  };
509 
510  template<typename ImageT, typename ThresholdT>
511  void findPeaks(PTR(detection::Footprint) foot, ImageT const& img, bool polarity, ThresholdT)
512  {
513  FindPeaksInFootprint<ImageT> peakFinder(img, polarity, foot->getPeaks());
514  peakFinder.apply(*foot, 1);
515 
516  // We use getInternal() here to get the vector of shared_ptr that Catalog uses internally,
517  // which causes the STL algorithm to copy pointers instead of PeakRecords (which is what
518  // it'd try to do if we passed Catalog's own iterators).
519  std::stable_sort(foot->getPeaks().getInternal().begin(), foot->getPeaks().getInternal().end(),
520  SortPeaks());
521 
522  if (foot->getPeaks().empty()) {
523  FindMaxInFootprint<ImageT> maxFinder(img, polarity);
524  maxFinder.apply(*foot);
525  maxFinder.addPeak(foot->getPeaks());
526  }
527  }
528 
529  // No need to search for peaks when processing a Mask
530  template<typename ImageT>
531  void findPeaks(PTR(detection::Footprint), ImageT const&, bool, ThresholdBitmask_traits)
532  {
533  ;
534  }
535 }
536 
537 /************************************************************************************************************/
538 /*
539  * Functions to determine if a pixel's in a Footprint
540  */
541 template<typename ImagePixelT, typename IterT>
542 static inline bool inFootprint(ImagePixelT pixVal, IterT,
543  bool polarity, double thresholdVal, ThresholdLevel_traits) {
544  return (polarity ? pixVal : -pixVal) >= thresholdVal;
545 }
546 
547 template<typename ImagePixelT, typename IterT>
548 static inline bool inFootprint(ImagePixelT pixVal, IterT var,
549  bool polarity, double thresholdVal, ThresholdPixelLevel_traits) {
550  return (polarity ? pixVal : -pixVal) >= thresholdVal*::sqrt(*var);
551 }
552 
553 template<typename ImagePixelT, typename IterT>
554 static inline bool inFootprint(ImagePixelT pixVal, IterT,
555  bool, double thresholdVal, ThresholdBitmask_traits) {
556  return (pixVal & static_cast<long>(thresholdVal));
557 }
558 
559 /*
560  * Advance the x_iterator to the variance image, when relevant (it may be NULL otherwise)
561  */
562 template<typename IterT>
563 static inline IterT
564 advancePtr(IterT varPtr, Threshold_traits) {
565  return varPtr;
566 }
567 
568 template<typename IterT>
569 static inline IterT
570 advancePtr(IterT varPtr, ThresholdPixelLevel_traits) {
571  return varPtr + 1;
572 }
573 
574 /*
575  * Here's the working routine for the FootprintSet constructors; see documentation
576  * of the constructors themselves
577  */
578 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT, typename ThresholdTraitT>
579 static void findFootprints(
580  typename detection::FootprintSet::FootprintList *_footprints, // Footprints
581  geom::Box2I const& _region, // BBox of pixels that are being searched
582  image::ImageBase<ImagePixelT> const &img, // Image to search for objects
583  image::Image<VariancePixelT> const *var, // img's variance
584  double const footprintThreshold, // threshold value for footprint
585  double const includeThresholdMultiplier, // threshold (relative to footprintThreshold) for inclusion
586  bool const polarity, // if false, search _below_ thresholdVal
587  int const npixMin, // minimum number of pixels in an object
588  bool const setPeaks // should I set the Peaks list?
589 )
590 {
591  int id; /* object ID */
592  int in_span; /* object ID of current IdSpan */
593  int nobj = 0; /* number of objects found */
594  int x0 = 0; /* unpacked from a IdSpan */
595 
596  typedef typename image::Image<ImagePixelT> ImageT;
597  double includeThreshold = footprintThreshold * includeThresholdMultiplier; // Threshold for inclusion
598 
599  int const row0 = img.getY0();
600  int const col0 = img.getX0();
601  int const height = img.getHeight();
602  int const width = img.getWidth();
603 /*
604  * Storage for arrays that identify objects by ID. We want to be able to
605  * refer to idp[-1] and idp[width], hence the (width + 2)
606  */
607  std::vector<int> id1(width + 2);
608  std::fill(id1.begin(), id1.end(), 0);
609  std::vector<int> id2(width + 2);
610  std::fill(id2.begin(), id2.end(), 0);
611  std::vector<int>::iterator idc = id1.begin() + 1; // object IDs in current/
612  std::vector<int>::iterator idp = id2.begin() + 1; // previous row
613 
614  std::vector<int> aliases; // aliases for initially disjoint parts of Footprints
615  aliases.reserve(1 + height/20); // initial size of aliases
616 
617  std::vector<IdSpan::Ptr> spans; // y:x0,x1 for objects
618  spans.reserve(aliases.capacity()); // initial size of spans
619 
620  aliases.push_back(0); // 0 --> 0
621 /*
622  * Go through image identifying objects
623  */
624  typedef typename image::Image<ImagePixelT>::x_iterator x_iterator;
625  typedef typename image::Image<VariancePixelT>::x_iterator x_var_iterator;
626 
627  in_span = 0; // not in a span
628  for (int y = 0; y != height; ++y) {
629  if (idc == id1.begin() + 1) {
630  idc = id2.begin() + 1;
631  idp = id1.begin() + 1;
632  } else {
633  idc = id1.begin() + 1;
634  idp = id2.begin() + 1;
635  }
636  std::fill_n(idc - 1, width + 2, 0);
637 
638  in_span = 0; /* not in a span */
639  bool good = (includeThresholdMultiplier == 1.0); /* Span exceeds the threshold? */
640 
641  x_iterator pixPtr = img.row_begin(y);
642  x_var_iterator varPtr = (var == NULL) ? NULL : var->row_begin(y);
643  for (int x = 0; x < width; ++x, ++pixPtr, varPtr = advancePtr(varPtr, ThresholdTraitT())) {
644  ImagePixelT const pixVal = *pixPtr;
645 
646  if (isBadPixel(pixVal) ||
647  !inFootprint(pixVal, varPtr, polarity, footprintThreshold, ThresholdTraitT())) {
648  if (in_span) {
649  IdSpan::Ptr sp(new IdSpan(in_span, y, x0, x - 1, good));
650  spans.push_back(sp);
651 
652  in_span = 0;
653  good = false;
654  }
655  } else { /* a pixel to fix */
656  if (idc[x - 1] != 0) {
657  id = idc[x - 1];
658  } else if (idp[x - 1] != 0) {
659  id = idp[x - 1];
660  } else if (idp[x] != 0) {
661  id = idp[x];
662  } else if (idp[x + 1] != 0) {
663  id = idp[x + 1];
664  } else {
665  id = ++nobj;
666  aliases.push_back(id);
667  }
668 
669  idc[x] = id;
670  if (!in_span) {
671  x0 = x;
672  in_span = id;
673  }
674 /*
675  * Do we need to merge ID numbers? If so, make suitable entries in aliases[]
676  */
677  if (idp[x + 1] != 0 && idp[x + 1] != id) {
678  aliases[resolve_alias(aliases, idp[x + 1])] = resolve_alias(aliases, id);
679 
680  idc[x] = id = idp[x + 1];
681  }
682 
683  if (!good && inFootprint(pixVal, varPtr, polarity, includeThreshold, ThresholdTraitT())) {
684  good = true;
685  }
686  }
687  }
688 
689  if (in_span) {
690  IdSpan::Ptr sp(new IdSpan(in_span, y, x0, width - 1, good));
691  spans.push_back(sp);
692  }
693  }
694 /*
695  * Resolve aliases; first alias chains, then the IDs in the spans
696  */
697  for (unsigned int i = 0; i < spans.size(); i++) {
698  spans[i]->id = resolve_alias(aliases, spans[i]->id);
699  }
700 /*
701  * Sort spans by ID, so we can sweep through them once
702  */
703  if (spans.size() > 0) {
704  std::sort(spans.begin(), spans.end(), IdSpanCompar());
705  }
706 /*
707  * Build Footprints from spans
708  */
709  unsigned int i0; // initial value of i
710  if (spans.size() > 0) {
711  id = spans[0]->id;
712  i0 = 0;
713  for (unsigned int i = 0; i <= spans.size(); i++) { // <= size to catch the last object
714  if (i == spans.size() || spans[i]->id != id) {
715  PTR(detection::Footprint) fp(new detection::Footprint(i - i0, _region));
716 
717  bool good = false; // Span includes pixel sufficient to include footprint in set?
718  for (; i0 < i; i0++) {
719  good |= spans[i0]->good;
720  fp->addSpan(spans[i0]->y + row0, spans[i0]->x0 + col0, spans[i0]->x1 + col0);
721  }
722 
723  if (good && !(fp->getNpix() < npixMin)) {
724  _footprints->push_back(fp);
725  }
726  }
727 
728  if (i < spans.size()) {
729  id = spans[i]->id;
730  }
731  }
732  }
733 /*
734  * Find all peaks within those Footprints
735  */
736  if (setPeaks) {
737  typedef detection::FootprintSet::FootprintList::iterator fiterator;
738  for (fiterator ptr = _footprints->begin(), end = _footprints->end(); ptr != end; ++ptr) {
739  findPeaks(*ptr, img, polarity, ThresholdTraitT());
740  }
741  }
742 }
743 
744 /************************************************************************************************************/
745 /*
746  * \brief Find a FootprintSet given an Image and a threshold
747  */
748 template<typename ImagePixelT>
750  image::Image<ImagePixelT> const &img,
751  Threshold const &threshold,
752  int const npixMin,
753  bool const setPeaks
754 ) : lsst::daf::base::Citizen(typeid(this)),
755  _footprints(new FootprintList()),
756  _region(img.getBBox())
757 {
758  typedef float VariancePixelT;
759 
760  findFootprints<ImagePixelT, afw::image::MaskPixel, VariancePixelT, ThresholdLevel_traits>(
761  _footprints.get(),
762  _region,
763  img,
764  NULL,
765  threshold.getValue(img), threshold.getIncludeMultiplier(), threshold.getPolarity(),
766  npixMin,
767  setPeaks
768  );
769 }
770 
771 // NOTE: not a template to appease swig (see note by instantiations at bottom)
772 
773 /*
774  * \brief Find a FootprintSet given a Mask and a threshold
775  */
776 template <typename MaskPixelT>
778  image::Mask<MaskPixelT> const &msk,
779  Threshold const &threshold,
780  int const npixMin
781 ) : lsst::daf::base::Citizen(typeid(this)),
782  _footprints(new FootprintList()),
783  _region(msk.getBBox())
784 {
785  switch (threshold.getType()) {
786  case Threshold::BITMASK:
787  findFootprints<MaskPixelT, MaskPixelT, float, ThresholdBitmask_traits>(
788  _footprints.get(), _region, msk, NULL, threshold.getValue(), threshold.getIncludeMultiplier(),
789  threshold.getPolarity(), npixMin, false);
790  break;
791 
792  case Threshold::VALUE:
793  findFootprints<MaskPixelT, MaskPixelT, float, ThresholdLevel_traits>(
794  _footprints.get(), _region, msk, NULL, threshold.getValue(), threshold.getIncludeMultiplier(),
795  threshold.getPolarity(), npixMin, false);
796  break;
797 
798  default:
799  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
800  "You must specify a numerical threshold value with a Mask");
801  }
802 }
803 
804 
817 template<typename ImagePixelT, typename MaskPixelT>
820  Threshold const &threshold,
821  std::string const &planeName,
822  int const npixMin,
823  bool const setPeaks
824 ) : lsst::daf::base::Citizen(typeid(this)),
825  _footprints(new FootprintList()),
826  _region(
827  geom::Point2I(maskedImg.getX0(), maskedImg.getY0()),
828  geom::Extent2I(maskedImg.getWidth(), maskedImg.getHeight())
829  )
830 {
831  typedef typename image::MaskedImage<ImagePixelT, MaskPixelT>::Variance::Pixel VariancePixelT;
832  // Find the Footprints
833  switch (threshold.getType()) {
835  findFootprints<ImagePixelT, MaskPixelT, VariancePixelT, ThresholdPixelLevel_traits>(
836  _footprints.get(),
837  _region,
838  *maskedImg.getImage(),
839  maskedImg.getVariance().get(),
840  threshold.getValue(maskedImg),
841  threshold.getIncludeMultiplier(),
842  threshold.getPolarity(),
843  npixMin,
844  setPeaks
845  );
846  break;
847  default:
848  findFootprints<ImagePixelT, MaskPixelT, VariancePixelT, ThresholdLevel_traits>(
849  _footprints.get(),
850  _region,
851  *maskedImg.getImage(),
852  maskedImg.getVariance().get(),
853  threshold.getValue(maskedImg),
854  threshold.getIncludeMultiplier(),
855  threshold.getPolarity(),
856  npixMin,
857  setPeaks
858  );
859  break;
860  }
861  // Set Mask if requested
862  if (planeName == "") {
863  return;
864  }
865  //
866  // Define the maskPlane
867  //
868  const typename image::Mask<MaskPixelT>::Ptr mask = maskedImg.getMask();
869  mask->addMaskPlane(planeName);
870 
871  MaskPixelT const bitPlane = mask->getPlaneBitMask(planeName);
872  //
873  // Set the bits where objects are detected
874  //
875  typedef image::Mask<MaskPixelT> MaskT;
876 
877  class MaskFootprint : public detection::FootprintFunctor<MaskT> {
878  public:
879  MaskFootprint(MaskT const& mimage,
880  MaskPixelT bit) : detection::FootprintFunctor<MaskT>(mimage), _bit(bit) {}
881 
882  void operator()(typename MaskT::xy_locator loc, int, int) {
883  *loc |= _bit;
884  }
885  private:
886  MaskPixelT _bit;
887  };
888 
889  MaskFootprint maskit(*maskedImg.getMask(), bitPlane);
890  for (FootprintList::const_iterator fiter = _footprints->begin();
891  fiter != _footprints->end(); ++fiter
892  ) {
893  Footprint::Ptr const foot = *fiter;
894 
895  maskit.apply(*foot);
896  }
897 }
898 
899 
900 /************************************************************************************************************/
901 namespace {
903  /*
904  * A data structure to hold the starting point for a search for pixels above threshold,
905  * used by pmFindFootprintAtPoint
906  *
907  * We don't want to find this span again --- it's already part of the footprint ---
908  * so we set appropriate mask bits
909  */
910  //
911  // An enum for what we should do with a Startspan
912  //
913  enum DIRECTION {DOWN = 0, // scan down from this span
914  UP, // scan up from this span
915  RESTART, // restart scanning from this span
916  DONE // this span is processed
917  };
918  //
919  // A Class that remembers how to [re-]start scanning the image for pixels
920  //
921  template<typename MaskPixelT>
922  class Startspan {
923  public:
924  typedef std::vector<boost::shared_ptr<Startspan> > Ptr;
925 
926  Startspan(detection::Span const *span, image::Mask<MaskPixelT> *mask, DIRECTION const dir);
927  ~Startspan() { delete _span; }
928 
929  bool getSpan() { return _span; }
930  bool Stop() { return _stop; }
931  DIRECTION getDirection() { return _direction; }
932 
933  static int detectedPlane; // The MaskPlane to use for detected pixels
934  static int stopPlane; // The MaskPlane to use for pixels that signal us to stop searching
935  private:
936  detection::Span::Ptr const _span; // The initial Span
937  DIRECTION _direction; // How to continue searching for further pixels
938  bool _stop; // should we stop searching?
939  };
940 
941  template<typename MaskPixelT>
942  Startspan<MaskPixelT>::Startspan(detection::Span const *span, // The span in question
943  image::Mask<MaskPixelT> *mask, // Pixels that we've already detected
944  DIRECTION const dir // Should we continue searching towards the top of the image?
945  ) :
946  _span(span),
947  _direction(dir),
948  _stop(false) {
949 
950  if (mask != NULL) { // remember that we've detected these pixels
951  mask->setMaskPlaneValues(detectedPlane, span->getX0(), span->getX1(), span->getY());
952 
953  int const y = span->getY() - mask->getY0();
954  for (int x = span->getX0() - mask->getX0(); x <= span->getX1() - mask->getX0(); x++) {
955  if (mask(x, y, stopPlane)) {
956  _stop = true;
957  break;
958  }
959  }
960  }
961  }
962 
963  template<typename ImagePixelT, typename MaskPixelT>
964  class StartspanSet {
965  public:
967  _image(image->getImage()),
968  _mask(image->getMask()) {}
969  ~StartspanSet() {}
970 
971  bool add(detection::Span *span, DIRECTION const dir, bool addToMask = true);
972  bool process(detection::Footprint *fp, // the footprint that we're building
973  detection::Threshold const &threshold, // Threshold
974  double const param = -1); // parameter that Threshold may need
975  private:
976  image::Image<ImagePixelT> const *_image; // the Image we're searching
977  image::Mask<MaskPixelT> *_mask; // the mask that tells us where we've got to
978  std::vector<typename Startspan<MaskPixelT>::Ptr> _spans; // list of Startspans
979  };
980 
981  //
982  // Add a new Startspan to a StartspansSet. Iff we see a stop bit, return true
983  //
984  template<typename ImagePixelT, typename MaskPixelT>
985  bool StartspanSet<ImagePixelT, MaskPixelT>::add(detection::Span *span, // the span in question
986  DIRECTION const dir, // the desired direction to search
987  bool addToMask) { // should I add the Span to the mask?
988  if (dir == RESTART) {
989  if (add(span, UP) || add(span, DOWN, false)) {
990  return true;
991  }
992  } else {
993  typename Startspan<MaskPixelT>::Ptr sspan(new Startspan<MaskPixelT>(span, dir));
994  if (sspan->stop()) { // we detected a stop bit
995  return true;
996  } else {
997  _spans.push_back(sspan);
998  }
999  }
1000 
1001  return false;
1002  }
1003 
1004  /************************************************************************************************/
1005  /*
1006  * Search the image for pixels above threshold, starting at a single Startspan.
1007  * We search the array looking for one to process; it'd be better to move the
1008  * ones that we're done with to the end, but it probably isn't worth it for
1009  * the anticipated uses of this routine.
1010  *
1011  * This is the guts of pmFindFootprintAtPoint
1012  */
1013  template<typename ImagePixelT, typename MaskPixelT>
1014  bool StartspanSet<ImagePixelT, MaskPixelT>::process(
1015  detection::Footprint *fp, // the footprint that we're building
1016  detection::Threshold const &threshold, // Threshold
1017  double const param // parameter that Threshold may need
1018  ) {
1019  int const row0 = _image->getY0();
1020  int const col0 = _image->getOffsetCols();
1021  int const height = _image->getHeight();
1022 
1023  /**********************************************************************************************/
1024 
1025  Startspan<MaskPixelT> *sspan = NULL;
1026  for (auto iter = _spans.begin(); iter != _spans.end(); iter++) {
1027  *sspan = *iter;
1028  if (sspan->getDirection() != DONE) {
1029  break;
1030  }
1031  if (sspan->Stop()) {
1032  break;
1033  }
1034  }
1035  if (sspan == NULL || sspan->getDirection() == DONE) { // no more Startspans to process
1036  return false;
1037  }
1038  if (sspan->Stop()) { // they don't want any more spans processed
1039  return false;
1040  }
1041  /*
1042  * Work
1043  */
1044  DIRECTION const dir = sspan->getDirection();
1045  /*
1046  * Set initial span to the startspan
1047  */
1048  int x0 = sspan->getSpan()->getX0() - col0;
1049  /*
1050  * Go through image identifying objects
1051  */
1052  int nx0 = -1; // new value of x0
1053  int const di = (dir == UP) ? 1 : -1; // how much i changes to get to the next row
1054  bool stop = false; // should I stop searching for spans?
1055 
1056  typedef typename image::Image<ImagePixelT>::pixel_accessor pixAccessT;
1057  double const thresholdVal = threshold.getValue(param);
1058  bool const polarity = threshold.getPolarity();
1059 
1060  for (int i = sspan->span->y -row0 + di; i < height && i >= 0; i += di) {
1061  pixAccessT imgRow = _image->origin().advance(0, i); // row pointer
1062  //maskPixAccessT maskRow = _mask->origin.advance(0, i); // masks's row pointer
1063  //
1064  // Search left from the pixel diagonally to the left of (i - di, x0). If there's
1065  // a connected span there it may need to grow up and/or down, so push it onto
1066  // the stack for later consideration
1067  //
1068  nx0 = -1;
1069  for (int j = x0 - 1; j >= -1; j--) {
1070  ImagePixelT pixVal = (j < 0) ? thresholdVal - 100 : (polarity ? imgRow[j] : -imgRow[j]);
1071  if (_mask(j, i, Startspan<MaskPixelT>::detectedPlane) || pixVal < threshold) {
1072  if (j < x0 - 1) { // we found some pixels above threshold
1073  nx0 = j + 1;
1074  }
1075  break;
1076  }
1077  }
1078 #if 0
1079  if (nx0 < 0) { // no span to the left
1080  nx1 = x0 - 1; // we're going to resume searching at nx1 + 1
1081  } else {
1082  //
1083  // Search right in leftmost span
1084  //
1085  //nx1 = 0; // make gcc happy
1086  for (int j = nx0 + 1; j <= width; j++) {
1087  ImagePixelT pixVal = (j >= width) ? threshold - 100 :
1088  (polarity ? (F32 ? imgRowF32[j] : imgRowS32[j]) :
1089  (F32 ? -imgRowF32[j] : -imgRowS32[j]));
1090  if ((maskRow[j] & DETECTED) || pixVal < threshold) {
1091  nx1 = j - 1;
1092  break;
1093  }
1094  }
1095 
1096  pmSpan const *sp = pmFootprintAddSpan(fp, i + row0, nx0 + col0, nx1 + col0);
1097 
1098  if (add_startspan(startspans, sp, mask, RESTART)) {
1099  stop = true;
1100  break;
1101  }
1102  }
1103  //
1104  // Now look for spans connected to the old span. The first of these we'll
1105  // simply process, but others will have to be deferred for later consideration.
1106  //
1107  // In fact, if the span overhangs to the right we'll have to defer the overhang
1108  // until later too, as it too can grow in both directions
1109  //
1110  // Note that column width exists virtually, and always ends the last span; this
1111  // is why we claim below that sx1 is always set
1112  //
1113  bool first = false; // is this the first new span detected?
1114  for (int j = nx1 + 1; j <= x1 + 1; j++) {
1115  ImagePixelT pixVal = (j >= width) ? threshold - 100 :
1116  (polarity ? (F32 ? imgRowF32[j] : imgRowS32[j]) : (F32 ? -imgRowF32[j] : -imgRowS32[j]));
1117  if (!(maskRow[j] & DETECTED) && pixVal >= threshold) {
1118  int sx0 = j++; // span that we're working on is sx0:sx1
1119  int sx1 = -1; // We know that if we got here, we'll also set sx1
1120  for (; j <= width; j++) {
1121  ImagePixelT pixVal = (j >= width) ? threshold - 100 :
1122  (polarity ? (F32 ? imgRowF32[j] : imgRowS32[j]) :
1123  (F32 ? -imgRowF32[j] : -imgRowS32[j]));
1124  if ((maskRow[j] & DETECTED) || pixVal < threshold) { // end of span
1125  sx1 = j;
1126  break;
1127  }
1128  }
1129  assert (sx1 >= 0);
1130 
1131  pmSpan const *sp;
1132  if (first) {
1133  if (sx1 <= x1) {
1134  sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);
1135  if (add_startspan(startspans, sp, mask, DONE)) {
1136  stop = true;
1137  break;
1138  }
1139  } else { // overhangs to right
1140  sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, x1 + col0);
1141  if (add_startspan(startspans, sp, mask, DONE)) {
1142  stop = true;
1143  break;
1144  }
1145  sp = pmFootprintAddSpan(fp, i + row0, x1 + 1 + col0, sx1 + col0 - 1);
1146  if (add_startspan(startspans, sp, mask, RESTART)) {
1147  stop = true;
1148  break;
1149  }
1150  }
1151  first = false;
1152  } else {
1153  sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);
1154  if (add_startspan(startspans, sp, mask, RESTART)) {
1155  stop = true;
1156  break;
1157  }
1158  }
1159  }
1160  }
1161 
1162  if (stop || first == false) { // we're done
1163  break;
1164  }
1165 
1166  x0 = nx0;
1167  x1 = nx1;
1168 #endif
1169  }
1170  /*
1171  * Cleanup
1172  */
1173 
1174  sspan->_direction = DONE;
1175  return stop ? false : true;
1176  }
1178 }
1179 #if 0
1180 
1181 
1182 /*
1183  * Go through an image, starting at (row, col) and assembling all the pixels
1184  * that are connected to that point (in a chess kings-move sort of way) into
1185  * a pmFootprint.
1186  *
1187  * This is much slower than pmFindFootprints if you want to find lots of
1188  * footprints, but if you only want a small region about a given point it
1189  * can be much faster
1190  *
1191  * N.b. The returned pmFootprint is not in "normal form"; that is the pmSpans
1192  * are not sorted by increasing y, x0, x1. If this matters to you, call
1193  * pmFootprintNormalize()
1194  */
1195 pmFootprint *
1196 pmFindFootprintAtPoint(psImage const *img, // image to search
1197  Threshold const &threshold, // Threshold
1198  psArray const *peaks, // array of peaks; finding one terminates search for footprint
1199  int row, int col) { // starting position (in img's parent's coordinate system)
1200  assert(img != NULL);
1201 
1202  bool F32 = false; // is this an F32 image?
1203  if (img->type.type == PS_TYPE_F32) {
1204  F32 = true;
1205  } else if (img->type.type == PS_TYPE_S32) {
1206  F32 = false;
1207  } else { // N.b. You can't trivially add more cases here; F32 is just a bool
1208  psError(PS_ERR_UNKNOWN, true, "Unsupported psImage type: %d", img->type.type);
1209  return NULL;
1210  }
1211  psF32 *imgRowF32 = NULL; // row pointer if F32
1212  psS32 *imgRowS32 = NULL; // " " " " !F32
1213 
1214  int const row0 = img->row0;
1215  int const col0 = img->col0;
1216  int const height = img->getHeight();
1217  int const width = img->getWidth();
1218 /*
1219  * Is point in image, and above threshold?
1220  */
1221  row -= row0;
1222  col -= col0;
1223  if (row < 0 || row >= height ||
1224  col < 0 || col >= width) {
1225  psError(PS_ERR_BAD_PARAMETER_VALUE, true,
1226  "row/col == (%d, %d) are out of bounds [%d--%d, %d--%d]",
1227  row + row0, col + col0, row0, row0 + height - 1, col0, col0 + width - 1);
1228  return NULL;
1229  }
1230 
1231  ImagePixelT pixVal = F32 ? img->data.F32[row][col] : img->data.S32[row][col];
1232  if (pixVal < threshold) {
1233  return pmFootprintAlloc(0, img);
1234  }
1235 
1236  pmFootprint *fp = pmFootprintAlloc(1 + img->getHeight()/10, img);
1237 /*
1238  * We need a mask for two purposes; to indicate which pixels are already detected,
1239  * and to store the "stop" pixels --- those that, once reached, should stop us
1240  * looking for the rest of the pmFootprint. These are generally set from peaks.
1241  */
1242  psImage *mask = psImageAlloc(width, height, PS_TYPE_MASK);
1243  P_PSIMAGE_SET_ROW0(mask, row0);
1244  P_PSIMAGE_SET_COL0(mask, col0);
1245  psImageInit(mask, INITIAL);
1246  //
1247  // Set stop bits from peaks list
1248  //
1249  assert (peaks == NULL || peaks->n == 0 || pmIsPeak(peaks->data[0]));
1250  if (peaks != NULL) {
1251  for (int i = 0; i < peaks->n; i++) {
1252  pmPeak *peak = peaks->data[i];
1253  mask->data.PS_TYPE_MASK_DATA[peak->y - mask->row0][peak->x - mask->col0] |= STOP;
1254  }
1255  }
1256 /*
1257  * Find starting span passing through (row, col)
1258  */
1259  psArray *startspans = psArrayAllocEmpty(1); // spans where we have to restart the search
1260 
1261  imgRowF32 = img->data.F32[row]; // only one of
1262  imgRowS32 = img->data.S32[row]; // these is valid!
1263  psMaskType *maskRow = mask->data.PS_TYPE_MASK_DATA[row];
1264  {
1265  int i;
1266  for (i = col; i >= 0; i--) {
1267  pixVal = F32 ? imgRowF32[i] : imgRowS32[i];
1268  if ((maskRow[i] & DETECTED) || pixVal < threshold) {
1269  break;
1270  }
1271  }
1272  int i0 = i;
1273  for (i = col; i < width; i++) {
1274  pixVal = F32 ? imgRowF32[i] : imgRowS32[i];
1275  if ((maskRow[i] & DETECTED) || pixVal < threshold) {
1276  break;
1277  }
1278  }
1279  int i1 = i;
1280  pmSpan const *sp = pmFootprintAddSpan(fp, row + row0, i0 + col0 + 1, i1 + col0 - 1);
1281 
1282  (void)add_startspan(startspans, sp, mask, RESTART);
1283  }
1284  /*
1285  * Now workout from those Startspans, searching for pixels above threshold
1286  */
1287  while (do_startspan(fp, img, mask, threshold, startspans)) continue;
1288  /*
1289  * Cleanup
1290  */
1291  psFree(mask);
1292  psFree(startspans); // restores the image pixel
1293 
1294  return fp; // pmFootprint really
1295 }
1296 #endif
1297 
1298 /************************************************************************************************************/
1303 ) :
1304  lsst::daf::base::Citizen(typeid(this)),
1305  _footprints(PTR(FootprintList)(new FootprintList)), _region(region) {
1306 }
1307 
1312  FootprintSet const &rhs
1313 ) :
1314  lsst::daf::base::Citizen(typeid(this)),
1315  _footprints(new FootprintList), _region(rhs._region)
1316 {
1317  _footprints->reserve(rhs._footprints->size());
1318  for (FootprintSet::FootprintList::const_iterator ptr = rhs._footprints->begin(),
1319  end = rhs._footprints->end(); ptr != end; ++ptr) {
1320  _footprints->push_back(PTR(Footprint)(new Footprint(**ptr)));
1321  }
1322 }
1323 
1327  FootprintSet tmp(rhs);
1328  swap(tmp); // See Meyers, Effective C++, Item 11
1329  return *this;
1330 }
1331 
1332 /************************************************************************************************************/
1337  detection::FootprintSet const& rhs,
1338  int tGrow,
1339  int rGrow,
1340  bool isotropic
1341 )
1342 {
1343  detection::FootprintControl const ctrl(true, isotropic);
1344  detection::FootprintSet fs = mergeFootprintSets(*this, tGrow, rhs, rGrow, ctrl);
1345  swap(fs); // Swap the new FootprintSet into place
1346 }
1347 
1351 //
1353  geom::Box2I const& region
1354 ) {
1355  _region = region;
1356 
1357  for (FootprintSet::FootprintList::iterator ptr = _footprints->begin(),
1358  end = _footprints->end(); ptr != end; ++ptr
1359  ) {
1360  (*ptr)->setRegion(region);
1361  }
1362 }
1363 
1364 /************************************************************************************************************/
1371  FootprintSet const &rhs,
1372  int r,
1373  bool isotropic
1374 )
1376  : lsst::daf::base::Citizen(typeid(this)), _footprints(new FootprintList), _region(rhs._region)
1377 {
1378  if (r == 0) {
1379  FootprintSet fs = rhs;
1380  swap(fs); // Swap the new FootprintSet into place
1381  return;
1382  } else if (r < 0) {
1383  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
1384  (boost::format("I cannot grow by negative numbers: %d") % r).str());
1385  }
1386 
1387  detection::FootprintControl const ctrl(true, isotropic);
1388  detection::FootprintSet fs = mergeFootprintSets(FootprintSet(rhs.getRegion()), 0, rhs, r, ctrl);
1389  swap(fs); // Swap the new FootprintSet into place
1390 }
1391 
1392 /************************************************************************************************************/
1393 
1395  int ngrow,
1396  detection::FootprintControl const& ctrl)
1397  : lsst::daf::base::Citizen(typeid(this)), _footprints(new FootprintList), _region(rhs._region)
1398 {
1399  if (ngrow == 0) {
1400  FootprintSet fs = rhs;
1401  swap(fs); // Swap the new FootprintSet into place
1402  return;
1403  } else if (ngrow < 0) {
1404  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
1405  str(boost::format("I cannot grow by negative numbers: %d") % ngrow));
1406  }
1407 
1408  detection::FootprintSet fs = mergeFootprintSets(FootprintSet(rhs.getRegion()), 0, rhs, ngrow, ctrl);
1409  swap(fs); // Swap the new FootprintSet into place
1410 }
1411 
1412 /************************************************************************************************************/
1419  FootprintSet const& fs1,
1420  FootprintSet const& fs2,
1421  bool const
1422  )
1423  : lsst::daf::base::Citizen(typeid(this)),
1424  _footprints(new FootprintList()),
1425  _region(fs1._region)
1426 {
1427  _region.include(fs2._region);
1428  throw LSST_EXCEPT(lsst::pex::exceptions::LogicError, "NOT IMPLEMENTED");
1429 }
1430 
1431 /************************************************************************************************************/
1438 detection::FootprintSet::insertIntoImage(
1439  bool const relativeIDs
1440 ) const {
1443  );
1444  *im = 0;
1445 
1447  for (FootprintList::const_iterator fiter = _footprints->begin();
1448  fiter != _footprints->end(); fiter++
1449  ) {
1450  Footprint::Ptr const foot = *fiter;
1451 
1452  if (relativeIDs) {
1453  id++;
1454  } else {
1455  id = foot->getId();
1456  }
1457 
1458  foot->insertIntoImage(*im.get(), id);
1459  }
1460 
1461  return im;
1462 }
1463 
1464 /************************************************************************************************************/
1468 template<typename ImagePixelT, typename MaskPixelT>
1469 void
1472  HeavyFootprintCtrl const *ctrl
1473 )
1474 {
1476 
1477  if (!ctrl) {
1478  ctrl = &ctrl_s;
1479  }
1480 
1481  for (FootprintList::iterator ptr = _footprints->begin(),
1482  end = _footprints->end(); ptr != end; ++ptr) {
1483  ptr->reset(new detection::HeavyFootprint<ImagePixelT, MaskPixelT>(**ptr, mimg, ctrl));
1484  }
1485 }
1486 
1489 ) const {
1490  for (FootprintList::const_iterator i = _footprints->begin(); i != _footprints->end(); ++i) {
1491  PTR(afw::table::SourceRecord) r = cat.addNew();
1492  r->setFootprint(*i);
1493  }
1494 }
1495 
1496 
1497 /************************************************************************************************************/
1498 //
1499 // Explicit instantiations
1500 //
1501 
1502 #ifndef DOXYGEN
1503 
1504 #define INSTANTIATE(PIXEL) \
1505  template detection::FootprintSet::FootprintSet( \
1506  image::Image<PIXEL> const &, Threshold const &, int const, bool const); \
1507  template detection::FootprintSet::FootprintSet( \
1508  image::MaskedImage<PIXEL,image::MaskPixel> const &, Threshold const &, \
1509  std::string const &, int const, bool const);\
1510  template void detection::FootprintSet::makeHeavy(image::MaskedImage<PIXEL,image::MaskPixel> const &, \
1511  HeavyFootprintCtrl const *)
1512 
1514  Threshold const &, int const);
1515 
1516 template void detection::FootprintSet::setMask(image::Mask<image::MaskPixel> *, std::string const &);
1517 template void detection::FootprintSet::setMask(PTR(image::Mask<image::MaskPixel>), std::string const &);
1518 
1519 INSTANTIATE(boost::uint16_t);
1520 INSTANTIATE(int);
1521 INSTANTIATE(float);
1522 INSTANTIATE(double);
1523 
1524 #endif // !DOXYGEN
int y
int iter
ImagePtr getImage(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s image.
Definition: MaskedImage.h:869
A Control Object for Footprints, controlling e.g. how they are grown.
Definition: FootprintCtrl.h:36
#define PTR(...)
Definition: base.h:41
x_iterator row_begin(int y) const
Definition: Image.h:319
for(FootprintList::const_iterator ptr=feet->begin(), end=feet->end();ptr!=end;++ptr)
Definition: saturated.cc:82
void swap(ImageBase< PixelT > &a, ImageBase< PixelT > &b)
Definition: Image.cc:291
A control object for HeavyFootprints.
Definition: FootprintCtrl.h:91
std::pair< bool, bool > isIsotropic() const
Return &lt;isSet, Value&gt; for isotropic grows.
Definition: FootprintCtrl.h:77
geom::Box2I const getRegion() const
Definition: FootprintSet.h:173
Internal & getInternal()
Return a reference to the internal vector-of-shared_ptr.
Definition: Catalog.h:658
#define CONST_PTR(...)
Definition: base.h:47
lsst::afw::detection::Footprint Footprint
Definition: Source.h:61
definition of the Trace messaging facilities
Use number of sigma given per-pixel s.d.
Definition: Threshold.h:52
Represent a set of pixels of an arbitrary shape and size, including values for those pixels; a HeavyF...
std::vector< Footprint::Ptr > FootprintList
The FootprintSet&#39;s set of Footprints.
Definition: FootprintSet.h:57
int getMinY() const
Definition: Box.h:125
A Threshold is used to pass a threshold value to detection algorithms.
Definition: Threshold.h:44
afw::table::CatalogT< PeakRecord > PeakCatalog
Definition: Peak.h:225
boost::shared_ptr< Footprint > Ptr
Definition: Footprint.h:67
Point< int, 2 > Point2I
Definition: PSF.h:39
iterator begin() const
Definition: Image.cc:325
The base class for all image classed (Image, Mask, MaskedImage, ...)
Definition: Image.h:115
static int addMaskPlane(const std::string &name)
Definition: Mask.cc:700
int isnan(T t)
Definition: ieee.h:110
int getIy() const
Convenience accessors for the keys in the minimal schema.
Definition: Peak.h:210
void setMaskPlaneValues(const int plane, const int x0, const int x1, const int y)
Set the bit specified by &quot;planeId&quot; for pixels (x0, y) ... (x1, y)
Definition: Mask.cc:1142
int const x0
Definition: saturated.cc:45
PixelT Pixel
A pixel in this ImageBase.
Definition: Image.h:133
boost::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
Definition: Catalog.h:469
boost::shared_ptr< Mask > Ptr
Definition: Mask.h:95
An integer coordinate rectangle.
Definition: Box.h:53
Vector param
A set of pixels in an Image, including those pixels&#39; actual values.
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
Represent a collections of footprints associated with image data.
VariancePtr getVariance(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s variance.
Definition: MaskedImage.h:890
void merge(FootprintSet const &rhs, int tGrow=0, int rGrow=0, bool isotropic=true)
boost::uint64_t FootprintIdPixel
Definition: FootprintSet.h:46
std::pair< bool, bool > isLeft() const
Definition: FootprintCtrl.h:65
std::pair< bool, bool > isCircular() const
Definition: FootprintCtrl.h:63
int getWidth() const
Return the number of columns in the image.
Definition: Image.h:237
boost::shared_ptr< Footprint > growFootprint(Footprint const &foot, int nGrow, bool isotropic=true)
Use (pixels &amp; (given mask))
Definition: Threshold.h:49
static MaskPixelT getPlaneBitMask(const std::vector< std::string > &names)
Return the bitmask corresponding to a vector of plane names OR&#39;d together.
Definition: Mask.cc:860
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
Definition: fwd.h:55
void setMask(image::Mask< MaskPixelT > *mask, std::string const &planeName)
Definition: FootprintSet.h:180
bool getPolarity() const
return Threshold&#39;s polarity
Definition: Threshold.h:88
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:93
float getPeakValue() const
Convenience accessors for the keys in the minimal schema.
Definition: Peak.h:219
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:77
if(width!=gim.getWidth()||height!=gim.getHeight()||x0!=gim.getX0()||y0!=gim.getY0())
Definition: saturated.cc:47
std::pair< bool, bool > isUp() const
Definition: FootprintCtrl.h:67
detection::FootprintSet::FootprintList FootprintList
Definition: saturated.cc:80
int getMinX() const
Definition: Box.h:124
A set of pixels in an Image.
Definition: Footprint.h:62
int getY0() const
Definition: Image.h:255
void setRegion(geom::Box2I const &region)
int x
double getIncludeMultiplier() const
return includeMultiplier
Definition: Threshold.h:94
int row
Definition: CR.cc:153
double getValue(const double param=-1) const
Definition: Threshold.cc:79
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
A set of Footprints, associated with a MaskedImage.
Definition: FootprintSet.h:53
std::pair< bool, bool > isRight() const
Definition: FootprintCtrl.h:66
#define INSTANTIATE(T)
Control Footprint-related algorithms.
MaskPtr getMask(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s mask.
Definition: MaskedImage.h:879
geom::Box2I _region
The corners of the MaskedImage that the detections live in.
Definition: FootprintSet.h:207
virtual void operator()(typename ImageT::xy_locator loc, int x, int y)=0
int id
Definition: CR.cc:151
int resolve_alias(const std::vector< int > &aliases, int id)
Definition: CR.cc:92
Compute Image Statistics.
afw::table::Key< double > b
void makeSources(afw::table::SourceCatalog &catalog) const
Add a new record corresponding to each footprint to a SourceCatalog.
Record class that contains measurements made on a single exposure.
Definition: Source.h:81
Implementation of the Class MaskedImage.
int getHeight() const
Return the number of rows in the image.
Definition: Image.h:239
std::pair< bool, bool > isDown() const
Definition: FootprintCtrl.h:68
FootprintSet(image::Image< ImagePixelT > const &img, Threshold const &threshold, int const npixMin=1, bool const setPeaks=true)
int const npixMin
Definition: saturated.cc:72
void makeHeavy(image::MaskedImage< ImagePixelT, MaskPixelT > const &mimg, HeavyFootprintCtrl const *ctrl=NULL)
Record class that represents a peak in a Footprint.
Definition: Peak.h:40
void swap(FootprintSet &rhs)
Definition: FootprintSet.h:130
ThresholdType getType() const
return type of threshold
Definition: Threshold.h:67
A functor class to allow users to process all the pixels in a Footprint.
FootprintSet & operator=(FootprintSet const &rhs)
Assignment operator.
ImageT val
Definition: CR.cc:154
A class to represent a 2-dimensional array of pixels.
Definition: PSF.h:43
boost::shared_ptr< FootprintList > _footprints
the Footprints of detected objects
Definition: FootprintSet.h:206
Extent< int, 2 > Extent2I
Definition: Extent.h:355
void insert(iterator pos, InputIterator first, InputIterator last, bool deep=false)
Insert an iterator range into the table.
Definition: Catalog.h:497
int col
Definition: CR.cc:152
void apply(Footprint const &foot, int const margin=0)
Apply operator() to each pixel in the Footprint.
boost::shared_ptr< FootprintList > getFootprints()
Definition: FootprintSet.h:146
int getIx() const
Convenience accessors for the keys in the minimal schema.
Definition: Peak.h:209
Include files required for standard LSST Exception handling.
int getX0() const
Definition: Image.h:247
size_type size() const
Return the number of elements in the catalog.
Definition: Catalog.h:401