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
Footprint.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2008-2015 LSST Corporation.
4  *
5  * This product includes software developed by the
6  * LSST Project (http://www.lsst.org/).
7  *
8  * This program is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the LSST License Statement and
19  * the GNU General Public License along with this program. If not,
20  * see <http://www.lsstcorp.org/LegalNotices/>.
21  */
22 
23 /*****************************************************************************/
28 #include <cassert>
29 #include <string>
30 #include <typeinfo>
31 #include <algorithm>
32 #include "boost/format.hpp"
33 #include "lsst/pex/logging/Trace.h"
34 #include "lsst/pex/exceptions.h"
35 #include "lsst/afw/image/Mask.h"
36 #include "lsst/afw/math/Kernel.h"
41 #include "lsst/afw/geom/Point.h"
46 #include "lsst/utils/ieee.h"
47 
48 namespace lsst {
49 namespace afw {
50 namespace detection {
51 
52 // anonymous namespace
53 namespace {
54 
55 /*
56  * Compare two Span%s by y, then x0, then x1
57  *
58  * A utility functor passed to sort; needed to dereference the boost::shared_ptrs.
59  */
60  struct compareSpanByYX :
61  public std::binary_function<Span::ConstPtr, Span::ConstPtr, bool> {
62  int operator()(Span::ConstPtr a, Span::ConstPtr b) {
63  return (*a) < (*b);
64  }
65  };
66 
68 geom::Point2D transformPoint(double x, double y,
69  image::Wcs const& source,
70  image::Wcs const& target){
71  return target.skyToPixel(*source.pixelToSky(x, y));
72 }
73 
74 
75 } //end namespace
76 
77 /*****************************************************************************/
79 int Footprint::id = 0;
80 
82  int nspan,
83  geom::Box2I const & region
84 ) : lsst::daf::base::Citizen(typeid(this)),
85  _fid(++id),
86  _area(0),
87  _bbox(geom::Box2I()),
88  _peaks(PeakTable::makeMinimalSchema()),
89  _region(region),
90  _normalized(true)
91 {
92  if (nspan < 0) {
93  throw LSST_EXCEPT(
94  lsst::pex::exceptions::InvalidParameterError,
95  str(boost::format("Number of spans requested is -ve: %d") % nspan));
96  }
97 }
98 
105  afw::table::Schema const & peakSchema,
106  int nspan,
107  geom::Box2I const & region
108 ) : lsst::daf::base::Citizen(typeid(this)),
109  _fid(++id),
110  _area(0),
111  _bbox(geom::Box2I()),
112  _peaks(peakSchema),
113  _region(region),
114  _normalized(true)
115 {
116  if (nspan < 0) {
117  throw LSST_EXCEPT(
118  lsst::pex::exceptions::InvalidParameterError,
119  str(boost::format("Number of spans requested is -ve: %d") % nspan));
120  }
121 }
126  geom::Box2I const& bbox,
127  geom::Box2I const& region
128 ) : lsst::daf::base::Citizen(typeid(this)),
129  _fid(++id),
130  _area(0),
131  _bbox(bbox),
132  _peaks(PeakTable::makeMinimalSchema()),
133  _region(region)
134 {
135  int const x0 = bbox.getMinX();
136  int const y0 = bbox.getMinY();
137  int const x1 = bbox.getMaxX();
138  int const y1 = bbox.getMaxY();
139 
140  for (int i = y0; i <= y1; ++i) {
141  addSpan(i, x0, x1);
142  }
143  _normalized=true;
144 }
145 
147  geom::Point2I const & center,
148  double const radius,
149  geom::BoxI const & region
150 ) : lsst::daf::base::Citizen(typeid(this)),
151  _fid(++id),
152  _area(0),
153  _bbox(geom::BoxI()),
154  _peaks(PeakTable::makeMinimalSchema()),
155  _region(region)
156 {
157  int const r2 = static_cast<int>(radius*radius + 0.5); // rounded radius^2
158  int const r = static_cast<int>(std::sqrt(static_cast<double>(r2))); // truncated radius; r*r <= r2
159 
160  for (int i = -r; i <= r; ++i) {
161  int hlen = static_cast<int>(std::sqrt(static_cast<double>(r2 - i*i)));
162  addSpan(center.getY() + i, center.getX() - hlen, center.getX() + hlen);
163  }
164  _normalized = true;
165 }
166 
168  geom::ellipses::Ellipse const & ellipse,
169  geom::Box2I const & region
170 ) : lsst::daf::base::Citizen(typeid(this)),
171  _fid(++id),
172  _area(0),
173  _bbox(geom::Box2I()),
174  _peaks(PeakTable::makeMinimalSchema()),
175  _region(region),
176  _normalized(true)
177 {
178  geom::ellipses::PixelRegion pr(ellipse);
179  for (
180  geom::ellipses::PixelRegion::Iterator spanIter = pr.begin(), end = pr.end();
181  spanIter != end;
182  ++spanIter
183  ) {
184  if (!spanIter->isEmpty()) {
185  addSpan(*spanIter);
186  }
187  }
188  _normalized = true;
189 }
190 
192  Footprint::SpanList const & spans,
193  geom::Box2I const & region
194 ) : lsst::daf::base::Citizen(typeid(this)),
195  _fid(++id),
196  _area(0),
197  _bbox(geom::Box2I()),
198  _peaks(PeakTable::makeMinimalSchema()),
199  _region(region),
200  _normalized(false)
201 {
202  _spans.reserve(spans.size());
203  for(SpanList::const_iterator i(spans.begin()); i != spans.end(); ++i) {
204  addSpan(**i);
205  }
206 }
207 
208 Footprint::Footprint(Footprint const & other)
209  : lsst::daf::base::Citizen(typeid(this)),
210  _fid(++id),
211  _bbox(other._bbox),
212  // peaks are deep-copied, but use the same Table as other
213  _peaks(other.getPeaks().getTable(), other.getPeaks().begin(), other.getPeaks().end(), true),
214  _region(other._region)
215 {
216  //deep copy spans
217  _spans.reserve(other._spans.size());
218  for(SpanList::const_iterator i(other._spans.begin());
219  i != other._spans.end(); ++i
220  ) {
221  addSpan(**i);
222  }
223  _area = other._area;
224  _normalized = other._normalized;
225 
226 }
227 
228 Footprint::~Footprint() {
229 }
230 
231 
232 PTR(PeakRecord) Footprint::addPeak(float fx, float fy, float value) {
233  PTR(PeakRecord) p = getPeaks().addNew();
234  p->setIx(fx);
235  p->setIy(fy);
236  p->setFx(fx);
237  p->setFy(fy);
238  p->setPeakValue(value);
239  return p;
240 }
241 
242 namespace {
243 
244 // comparison function to sort peaks from most positive to most negative.
245 struct SortPeaks {
246 
247  SortPeaks(afw::table::Key<float> const & key) : _key(key) {}
248 
249  bool operator()(detection::PeakRecord const & a, detection::PeakRecord const & b) const {
250  return a.get(_key) > b.get(_key);
251  }
252 
253 private:
254  afw::table::Key<float> _key;
255 };
256 
257 } // anonymous
258 
259 void Footprint::sortPeaks(afw::table::Key<float> const & key) {
260  getPeaks().sort(SortPeaks(key.isValid() ? key : PeakTable::getPeakValueKey()));
261 }
262 
266 bool Footprint::contains(
267  lsst::afw::geom::Point2I const& pix
268 ) const
269 {
270  if (_bbox.contains(pix)) {
271  for (Footprint::SpanList::const_iterator siter = _spans.begin(); siter != _spans.end(); ++siter) {
272  if ((*siter)->contains(pix.getX(), pix.getY())) {
273  return true;
274  }
275  }
276  }
277 
278  return false;
279 }
280 
281 namespace {
282 
284 struct ClipSpansPredicate : public std::unary_function<PTR(Span) const&, bool> {
285  geom::Box2I const& bbox;
286 
287  ClipSpansPredicate(geom::Box2I const& _bbox) : bbox(_bbox) {}
288 
289  bool operator()(PTR(Span) const& spanPtr) const {
290  Span const& span = *spanPtr;
291  int const y = span.getY();
292  return (y < bbox.getMinY() || y > bbox.getMaxY() ||
293  span.getX0() > bbox.getMaxX() || span.getX1() < bbox.getMinX());
294  }
295 };
296 
298 struct ClipPeaksPredicate : public std::unary_function<PeakRecord const&, bool> {
299  geom::Box2I const& bbox;
300 
301  ClipPeaksPredicate(geom::Box2I const& _bbox) : bbox(_bbox) {}
302 
303  bool operator()(PTR(PeakRecord) const& peak) const {
304  return !bbox.contains(geom::Point2I(peak->getIx(), peak->getIy()));
305  }
306 };
307 } // anonymous namespace
308 
309 void Footprint::clipTo(geom::Box2I const& bbox) {
310  _spans.erase(std::remove_if(_spans.begin(), _spans.end(), ClipSpansPredicate(bbox)), _spans.end());
311  for (SpanList::const_iterator ss = _spans.begin(); ss != _spans.end(); ++ss) {
312  Span& span = **ss;
313  span.getX0() = std::max(span.getX0(), bbox.getMinX());
314  span.getX1() = std::min(span.getX1(), bbox.getMaxX());
315  }
316 
317  // Remove peaks not in the new bbox
318  _peaks.getInternal().erase(std::remove_if(_peaks.getInternal().begin(), _peaks.getInternal().end(),
319  ClipPeaksPredicate(bbox)), _peaks.getInternal().end());
320 
321  if (_spans.empty()) {
322  _bbox = geom::Box2I();
323  _area = 0;
324  _normalized = true;
325  } else {
326  _normalized = false;
327  normalize();
328  }
329 }
330 
331 void Footprint::normalize() {
332  if (_normalized) {
333  return;
334  }
335  assert(!_spans.empty());
336  //
337  // Check that the spans are sorted, and (more importantly) that each pixel appears
338  // in only one span
339  //
340  sort(_spans.begin(), _spans.end(), compareSpanByYX());
341 
342  Footprint::SpanList::iterator ptr = _spans.begin(), end = _spans.end();
343  Span *lspan = ptr->get(); // Left span
344  int y = lspan->_y;
345  int x1 = lspan->_x1;
346  _area = lspan->getWidth();
347  int minX = lspan->_x0, minY=y, maxX=x1;
348 
349  ++ptr;
350 
351  for (; ptr != end; ++ptr) {
352  Span *rspan = ptr->get(); // Right span
353  if (rspan->_y == y) {
354  if (rspan->_x0 <= x1 + 1) { // Spans overlap or touch
355  if (rspan->_x1 > x1) { // right span extends left span
356  //update area
357  _area += rspan->_x1 - x1;
358  //update end of current span
359  x1 = lspan->_x1 = rspan->_x1;
360  //update bounds
361  if(x1 > maxX) maxX = x1;
362  }
363 
364  ptr = _spans.erase(ptr);
365  end = _spans.end(); // delete the right span
366  if (ptr == end) {
367  break;
368  }
369 
370  --ptr;
371  continue;
372  }
373  else{
374  _area += rspan->getWidth();
375  if(rspan->_x1 > maxX) maxX = rspan->_x1;
376  }
377  } else {
378  _area += rspan->getWidth();
379  }
380 
381  y = rspan->_y;
382  x1 = rspan->_x1;
383 
384  lspan = rspan;
385  if(lspan->_x0 < minX) minX = lspan->_x0;
386  if(x1 > maxX) maxX = x1;
387  }
388  _bbox = geom::Box2I(geom::Point2I(minX, minY), geom::Point2I(maxX, y));
389 
390  _normalized = true;
391 }
392 
393 Span const& Footprint::addSpan(
394  int const y,
395  int const x0,
396  int const x1
397 ) {
398  if (x1 < x0) {
399  return addSpan(y, x1, x0);
400  }
401 
402  PTR(Span) sp(new Span(y, x0, x1));
403  _spans.push_back(sp);
404 
405  _area += sp->getWidth();
406  _normalized = false;
407 
408  _bbox.include(geom::Point2I(x0, y));
409  _bbox.include(geom::Point2I(x1, y));
410 
411  return *sp.get();
412 }
413 
414 const Span& Footprint::addSpan(
415  Span const& span
416 ) {
417  return addSpan(span._y, span._x0, span._x1);
418 }
419 
420 const Span& Footprint::addSpan(
421  Span const& span,
422  int dx,
423  int dy
424 ) {
425  return addSpan(span._y + dy, span._x0 + dx, span._x1 + dx);
426 }
427 
428 const Span& Footprint::addSpanInSeries(
429  int const y,
430  int const x0,
431  int const x1
432 ) {
433  if (x1 < x0) {
434  return addSpanInSeries(y, x1, x0);
435  }
436  if (_spans.size() == 0) {
437  const Span& s = addSpan(y, x0, x1);
438  _normalized = true;
439  return s;
440  }
441  // merge contiguous spans
442  PTR(Span) lastspan = _spans.back();
443  if ((y == lastspan->getY()) &&
444  (x0 == (lastspan->getX1() + 1))) {
445  // contiguous.
446  lastspan->_x1 = x1;
447  _area += (1 + x1 - x0);
448  _bbox.include(geom::Point2I(x1,y));
449  return *lastspan;
450  }
451  if (!((y > lastspan->getY()) ||
452  (x0 > (lastspan->getX1() + 1)))) {
453  throw LSST_EXCEPT(
454  lsst::pex::exceptions::InvalidParameterError,
455  str(boost::format("addSpanInSeries: new span %i,[%i,%i] is NOT in series after last span "
456  "%i,[%i,%i]") %
457  y % x0 % x1 % lastspan->getY() % lastspan->getX0() % lastspan->getX1()));
458  }
459  const Span& s = addSpan(y, x0, x1);
460  _normalized = true;
461  return s;
462 }
463 
464 void Footprint::shift(
465  int dx,
466  int dy
467 ) {
468  for (SpanList::iterator i = _spans.begin(); i != _spans.end(); ++i){
469  PTR(Span) span = *i;
470 
471  span->_y += dy;
472  span->_x0 += dx;
473  span->_x1 += dx;
474  }
475 
476  _bbox.shift(geom::Extent2I(dx, dy));
477 }
478 
479 geom::Point2D
480 Footprint::getCentroid() const
481 {
482  int n = 0;
483  double xc = 0, yc = 0;
484  for (Footprint::SpanList::const_iterator siter = _spans.begin(); siter != _spans.end(); ++siter) {
485  CONST_PTR(Span) span = *siter;
486  int const y = span->getY();
487  int const x0 = span->getX0();
488  int const x1 = span->getX1();
489  int const npix = x1 - x0 + 1;
490 
491  n += npix;
492  xc += npix*0.5*(x1 + x0);
493  yc += npix*y;
494  }
495  assert(n == _area);
496 
497  return geom::Point2D(xc/_area, yc/_area);
498 }
499 
500 geom::ellipses::Quadrupole
501 Footprint::getShape() const
502 {
503  geom::Point2D cen = getCentroid();
504  double const xc = cen.getX();
505  double const yc = cen.getY();
506 
507  double sumxx = 0, sumxy = 0, sumyy = 0;
508  for (Footprint::SpanList::const_iterator siter = _spans.begin(); siter != _spans.end(); ++siter) {
509  CONST_PTR(Span) span = *siter;
510  int const y = span->getY();
511  int const x0 = span->getX0();
512  int const x1 = span->getX1();
513  int const npix = x1 - x0 + 1;
514 
515  for (int x = x0; x <= x1; ++x) {
516  sumxx += (x - xc)*(x - xc);
517  }
518  sumxy += npix*(0.5*(x1 + x0) - xc)*(y - yc);
519  sumyy += npix*(y - yc)*(y - yc);
520  }
521 
522  return geom::ellipses::Quadrupole(sumxx/_area, sumyy/_area, sumxy/_area);
523 }
524 
525 namespace {
526  /*
527  * Set the pixels in idImage which are in Footprint by adding or
528  * replacing the specified value to the Image
529  *
530  * The ids that are overwritten are returned for the callers deliction
531  */
532  template<bool overwriteId, typename PixelT>
533  void
534  doInsertIntoImage(geom::Box2I const& _region, // unpacked from Footprint
535  Footprint::SpanList const& _spans, // unpacked from Footprint
536  image::Image<PixelT>& idImage, // Image to contain the footprint
537  boost::uint64_t const id, // Add/replace id to idImage for pixels in Footprint
538  geom::Box2I const& region, // Footprint's region (default: getRegion())
539  long const mask=0x0, // Don't overwrite bits in this mask
540  std::set<boost::uint64_t> *oldIds=NULL // if non-NULL, set the IDs that were overwritten
541  )
542  {
543  int width, height, x0, y0;
544  if(!region.isEmpty()) {
545  height = region.getHeight();
546  width = region.getWidth();
547  x0 = region.getMinX();
548  y0 = region.getMinY();
549  } else {
550  height = _region.getHeight();
551  width = _region.getWidth();
552  x0 = _region.getMinX();
553  y0 = _region.getMinY();
554  }
555 
556  if (width != idImage.getWidth() || height != idImage.getHeight()) {
557  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
558  str(boost::format("Image of size (%dx%d) doesn't match "
559  "Footprint's host Image of size (%dx%d)") %
560  idImage.getWidth() % idImage.getHeight() % width % height));
561  }
562 
563  if (id & mask) {
564  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
565  str(boost::format("Id 0x%x sets bits in the protected mask 0x%x") % id % mask));
566  }
567 
568  typename std::set<boost::uint64_t>::const_iterator pos; // hint on where to insert into oldIds
569  if (oldIds) {
570  pos = oldIds->begin();
571  }
572  for (Footprint::SpanList::const_iterator spi = _spans.begin(); spi != _spans.end(); ++spi) {
573  CONST_PTR(Span) span = *spi;
574 
575  int const sy0 = span->getY() - y0;
576  if (sy0 < 0 || sy0 >= height) {
577  continue;
578  }
579 
580  int sx0 = span->getX0() - x0;
581  if (sx0 < 0) {
582  sx0 = 0;
583  }
584  int sx1 = span->getX1() - x0;
585  int const swidth = (sx1 >= width) ? width - sx0 : sx1 - sx0 + 1;
586 
587  for (typename image::Image<PixelT>::x_iterator ptr = idImage.x_at(sx0, sy0),
588  end = ptr + swidth; ptr != end; ++ptr) {
589  if (overwriteId) {
590  long val = *ptr & ~mask;
591  if (val != 0 and oldIds != NULL) {
592  pos = oldIds->insert(pos, val); // update our hint, pos
593  }
594  *ptr = (*ptr & mask) + id;
595  } else {
596  *ptr += id;
597  }
598  }
599  }
600  }
601 }
602 
603 template<typename PixelT>
604 void
605 Footprint::clipToNonzero(typename image::Image<PixelT> const& img) {
606  typedef lsst::afw::image::Image<PixelT> ImageT;
607  int const ix0 = img.getX0();
608  int const iy0 = img.getY0();
609  PixelT const zero = 0;
610 
611  normalize(); // allows us to produce a normalized output
612  SpanList old;
613  std::swap(_spans, old);
614  _spans.reserve(old.size());
615  _area = 0;
616  _bbox = geom::Box2I();
617  for (SpanList::iterator s = old.begin(); s != old.end(); ++s) {
618  int const y = (*s)->getY();
619  int const x0 = (*s)->getX0();
620  int const x1 = (*s)->getX1();
621  typename ImageT::x_iterator img_it = img.row_begin(y - iy0) + (x0 - ix0);
622  int leftx, rightx;
623  // find zero pixels on the left...
624  for (leftx = x0; leftx <= x1; ++leftx, ++img_it) {
625  if (*img_it != zero) {
626  break;
627  }
628  }
629  if (leftx > x1) {
630  // whole span is zero; drop it.
631  continue;
632  }
633  // find zero pixels on the right...
634  img_it = img.row_begin(y - iy0) + (x1 - ix0);
635  for (rightx = x1; rightx >= leftx; --rightx, --img_it) {
636  if (*img_it != zero) {
637  break;
638  }
639  }
640  addSpanInSeries(y, leftx, rightx);
641  }
642  normalize();
643 }
644 
645 template<typename PixelT>
646 void
647 Footprint::insertIntoImage(
648  typename image::Image<PixelT>& idImage,
649  boost::uint64_t const id,
650  geom::Box2I const& region
651 ) const
652 {
653  static_cast<void>(doInsertIntoImage<false>(_region, _spans, idImage, id, region));
654 }
655 
656 template<typename PixelT>
657 void
658 Footprint::insertIntoImage(
659  image::Image<PixelT>& idImage,
660  boost::uint64_t const id,
661  bool overwriteId,
662  long const mask,
663  std::set<boost::uint64_t> *oldIds,
664  geom::Box2I const& region
665 ) const
666 {
667  if (id > std::size_t(std::numeric_limits<PixelT>::max())) {
668  throw LSST_EXCEPT(
669  lsst::pex::exceptions::OutOfRangeError,
670  "id out of range for image type"
671  );
672  }
673  if (overwriteId) {
674  doInsertIntoImage<true>(_region, _spans, idImage, id, region, mask, oldIds);
675  } else {
676  doInsertIntoImage<false>(_region, _spans, idImage, id, region, mask, oldIds);
677  }
678 }
679 
680 void Footprint::include(std::vector<PTR(Footprint)> const & others, bool ignoreSelf) {
681  if (others.empty()) return;
682  geom::Box2I bbox;
683  if (!ignoreSelf) {
684  bbox.include(getBBox());
685  } else {
686  _spans.clear();
687  }
688  for (std::vector<PTR(Footprint)>::const_iterator i = others.begin(); i != others.end(); ++i) {
689  bbox.include((**i).getBBox());
690  }
691  boost::uint16_t bits = 0x1;
692  image::Mask<boost::uint16_t> mask(bbox);
693  if (!ignoreSelf) {
694  setMaskFromFootprint(&mask, *this, bits);
695  }
696  for (std::vector<PTR(Footprint)>::const_iterator i = others.begin(); i != others.end(); ++i) {
697  setMaskFromFootprint(&mask, **i, bits);
698  }
699  FootprintSet fpSet(mask, Threshold(bits, Threshold::BITMASK));
700  if (fpSet.getFootprints()->empty()) {
701  _spans.clear();
702  } else if (fpSet.getFootprints()->size() == 1u) {
703  _spans.swap(fpSet.getFootprints()->front()->getSpans());
704  } else {
705  _spans.clear();
706  for (std::vector<PTR(Footprint)>::const_iterator i = fpSet.getFootprints()->begin();
707  i != fpSet.getFootprints()->end(); ++i) {
708  _spans.insert(_spans.end(), (**i).getSpans().begin(), (**i).getSpans().end());
709  }
710  }
711  _normalized = false;
712  normalize();
713 }
714 
715 // Factory class used for table-based persistence; invoked via registry in afw::table::io
716 class FootprintFactory : public table::io::PersistableFactory {
717 public:
718 
719  virtual PTR(table::io::Persistable)
720  read(InputArchive const & archive, CatalogVector const & catalogs) const {
721  LSST_ARCHIVE_ASSERT(catalogs.size() == 2u);
722  PTR(Footprint) result = boost::make_shared<Footprint>();
723  result->readSpans(catalogs.front());
724  result->readPeaks(catalogs.back());
725  return result;
726  }
727 
728  explicit FootprintFactory(std::string const & name) : table::io::PersistableFactory(name) {}
729 
730 };
731 
732 namespace {
733 
734 // Singleton helper class that manages the schema and keys for persisting a Footprint
735 class FootprintPersistenceHelper : private boost::noncopyable {
736 public:
737  table::Schema spanSchema;
738  table::Key<int> spanY;
739  table::Key<int> spanX0;
740  table::Key<int> spanX1;
741 
742  static FootprintPersistenceHelper const & get() {
743  static FootprintPersistenceHelper instance;
744  return instance;
745  }
746 
747 private:
748  FootprintPersistenceHelper() :
749  spanSchema(),
750  spanY(spanSchema.addField<int>("y", "row position of span", "pixels")),
751  spanX0(spanSchema.addField<int>("x0", "first column of span (inclusive)", "pixels")),
752  spanX1(spanSchema.addField<int>("x1", "first column of span (inclusive)", "pixels"))
753  {
754  spanSchema.getCitizen().markPersistent();
755  }
756 };
757 
758 std::string getFootprintPersistenceName() { return "Footprint"; }
759 
760 // Insert the factory into the registry (instantiating an instance is sufficient, because
761 // the code that does the work is in the base class ctor)
762 FootprintFactory registration(getFootprintPersistenceName());
763 
764 } // anonymous
765 
766 std::string Footprint::getPersistenceName() const { return getFootprintPersistenceName(); }
767 
768 std::string Footprint::getPythonModule() const { return "lsst.afw.detection"; }
769 
770 void Footprint::write(OutputArchiveHandle & handle) const {
771  FootprintPersistenceHelper const & keys = FootprintPersistenceHelper::get();
772  afw::table::BaseCatalog spanCat = handle.makeCatalog(keys.spanSchema);
773  spanCat.reserve(_spans.size());
774  for (SpanList::const_iterator i = _spans.begin(); i != _spans.end(); ++i) {
775  PTR(afw::table::BaseRecord) record = spanCat.addNew();
776  record->set(keys.spanY, (**i).getY());
777  record->set(keys.spanX0, (**i).getX0());
778  record->set(keys.spanX1, (**i).getX1());
779  }
780  handle.saveCatalog(spanCat);
781  afw::table::BaseCatalog peakCat = handle.makeCatalog(_peaks.getSchema());
782  peakCat.insert(peakCat.end(), _peaks.begin(), _peaks.end(), true);
783  handle.saveCatalog(peakCat);
784 }
785 
786 void Footprint::readSpans(afw::table::BaseCatalog const & spanCat) {
787  FootprintPersistenceHelper const & keys = FootprintPersistenceHelper::get();
788  for (afw::table::BaseCatalog::const_iterator i = spanCat.begin(); i != spanCat.end(); ++i) {
789  addSpan(i->get(keys.spanY), i->get(keys.spanX0), i->get(keys.spanX1));
790  }
791 }
792 
793 void Footprint::readPeaks(afw::table::BaseCatalog const & peakCat) {
794  if (!peakCat.getSchema().contains(PeakTable::makeMinimalSchema())) {
795  // need to handle an older form of Peak persistence for backwards compatibility
796  afw::table::SchemaMapper mapper(peakCat.getSchema());
797  mapper.addMinimalSchema(PeakTable::makeMinimalSchema());
798  afw::table::Key<float> oldX = peakCat.getSchema()["x"];
799  afw::table::Key<float> oldY = peakCat.getSchema()["y"];
800  afw::table::Key<float> oldPeakValue = peakCat.getSchema()["value"];
801  mapper.addMapping(oldX, "f.x");
802  mapper.addMapping(oldY, "f.y");
803  mapper.addMapping(oldPeakValue, "peakValue");
804  _peaks = PeakCatalog(mapper.getOutputSchema());
805  _peaks.reserve(peakCat.size());
806  for (afw::table::BaseCatalog::const_iterator i = peakCat.begin(); i != peakCat.end(); ++i) {
807  PTR(PeakRecord) newPeak = _peaks.addNew();
808  newPeak->assign(*i, mapper);
809  newPeak->setIx(int(newPeak->getFx()));
810  newPeak->setIy(int(newPeak->getFy()));
811  }
812  return;
813  }
814  _peaks = PeakCatalog(peakCat.getSchema());
815  _peaks.reserve(peakCat.size());
816  for (afw::table::BaseCatalog::const_iterator i = peakCat.begin(); i != peakCat.end(); ++i) {
817  _peaks.addNew()->assign(*i);
818  }
819 }
820 
824 Footprint & Footprint::operator=(Footprint & other) {
825  _region = other._region;
826 
827  //deep copy spans
828  _spans = SpanList();
829  _spans.reserve(other._spans.size());
830  for(SpanList::const_iterator i(other._spans.begin());
831  i != other._spans.end(); ++i
832  ) {
833  addSpan(**i);
834  }
835  _area = other._area;
836  _normalized = other._normalized;
837  _bbox = other._bbox;
838 
839  //deep copy peaks
840  _peaks = PeakCatalog(other.getPeaks().getTable(), other.getPeaks().begin(), other.getPeaks().end(), true);
841  return *this;
842 }
843 
844 template<typename MaskT>
845 void Footprint::intersectMask(
846  lsst::afw::image::Mask<MaskT> const & mask,
847  MaskT const bitmask
848 ) {
849  geom::Box2I maskBBox = mask.getBBox();
850 
851  //this operation makes no sense on non-normalized footprints.
852  //make sure this is normalized
853  normalize();
854 
855  SpanList::iterator s(_spans.begin());
856  while((*s)->getY() < maskBBox.getMinY() && s != _spans.end()){
857  ++s;
858  }
859 
860 
861  int x0, x1, y;
862  SpanList maskedSpans;
863  int maskedArea=0;
864  for( ; s != _spans.end(); ++s) {
865  y = (*s)->getY();
866 
867  if (y > maskBBox.getMaxY())
868  break;
869 
870  x0 = (*s)->getX0();
871  x1 = (*s)->getX1();
872 
873  if(x1 < maskBBox.getMinX() || x0 > maskBBox.getMaxX()) {
874  //span is entirely outside the image mask. cannot be used
875  continue;
876  }
877 
878  //clip the span to be within the mask
879  if(x0 < maskBBox.getMinX()) x0 = maskBBox.getMinX();
880  if(x1 > maskBBox.getMaxX()) x1 = maskBBox.getMaxX();
881 
882  //Image iterators are always specified with respect to (0,0)
883  //regardless what the image::XY0 is set to.
884  typename image::Mask<MaskT>::const_x_iterator mIter = mask.x_at(
885  x0 - maskBBox.getMinX(), y - maskBBox.getMinY()
886  );
887 
888  //loop over all span locations, slicing the span at maskedPixels
889  for(int x = x0; x <= x1; ++x, ++mIter) {
890  if((*mIter & bitmask) != 0) {
891  //masked pixel found within span
892  if (x > x0) {
893  //add beginning of span to the output
894  //the fixed span contains all the unmasked pixels up to,
895  //but not including this masked pixel
896  PTR(Span) maskedSpan(new Span(y, x0, x- 1));
897  maskedSpans.push_back(maskedSpan);
898  maskedArea += maskedSpan->getWidth();
899  }
900  //set the next Span to start after this pixel
901  x0 = x + 1;
902  }
903  }
904 
905  //add last section of span
906  if(x0 <= x1) {
907  PTR(Span) maskedSpan(new Span(y, x0, x1));
908  maskedSpans.push_back(maskedSpan);
909  maskedArea += maskedSpan->getWidth();
910  }
911  }
912  _area = maskedArea;
913  _spans = maskedSpans;
914  _bbox.clip(maskBBox);
915 }
916 
917 
918 PTR(Footprint) Footprint::transform(
919  image::Wcs const& source,
920  image::Wcs const& target,
921  geom::Box2I const& region,
922  bool doClip
923 ) const {
924  // Transform the original bounding box
925  geom::Box2I const& fpBox = getBBox(); // Original bounding box
926  geom::Box2D tBoxD;
927  // If slow, could consider linearising the WCSes and combining the
928  // linear versions to a single transform, and then using that to
929  // transform all the points.
930  tBoxD.include(transformPoint(fpBox.getMinX(), fpBox.getMinY(), source, target));
931  tBoxD.include(transformPoint(fpBox.getMinX(), fpBox.getMaxY(), source, target));
932  tBoxD.include(transformPoint(fpBox.getMaxX(), fpBox.getMinY(), source, target));
933  tBoxD.include(transformPoint(fpBox.getMaxX(), fpBox.getMaxY(), source, target));
934  geom::Box2I tBoxI(tBoxD);
935 
936  // enumerate points in the new bbox that, when reverse-transformed, are within the given footprint.
937  PTR(Footprint) fpNew = boost::make_shared<Footprint>(getPeaks().getSchema(), 0, region);
938 
939  for (int y = tBoxI.getBeginY(); y < tBoxI.getEndY(); ++y) {
940  bool inSpan = false; // Are we in a span?
941  int start = -1; // Start of span
942 
943  for (int x = tBoxI.getBeginX(); x < tBoxI.getEndX(); ++x) {
944  geom::Point2D p = transformPoint(x, y, target, source);
945  int const xSource = std::floor(0.5 + p.getX());
946  int const ySource = std::floor(0.5 + p.getY());
947 
948  if (contains(geom::Point2I(xSource, ySource))) {
949  if (!inSpan) {
950  inSpan = true;
951  start = x;
952  }
953  } else if (inSpan) {
954  inSpan = false;
955  fpNew->addSpan(y, start, x - 1);
956  }
957  }
958  if (inSpan) {
959  fpNew->addSpan(y, start, tBoxI.getMaxX());
960  }
961  }
962 
963  // Copy over peaks to new Footprint
964  for (
965  PeakCatalog::const_iterator iter = this->getPeaks().begin();
966  iter != this->getPeaks().end();
967  ++iter
968  ) {
969  geom::Point2D tp = transformPoint(iter->getFx(), iter->getFy(), source, target);
970  fpNew->addPeak(tp.getX(), tp.getY(), iter->getPeakValue());
971  }
972 
973  if (doClip) {
974  fpNew->clipTo(region);
975  }
976  return fpNew;
977 }
978 
979 PTR(Footprint) Footprint::findEdgePixels() const
980 {
981  if (!_normalized) {
982  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "Footprint isn't normalized");
983  }
984  int const width = getBBox().getWidth(), height = getBBox().getHeight();
985  if (height <= 2 || _spans.size() <= 2) {
986  // Everything is on the edge
987  return boost::make_shared<Footprint>(*this);
988  }
989 
990  // Get a list of pixels (in the form of a Footprint) that are on the edge horizontally
991  // or have nothing above or below them.
992  PTR(Footprint) edges = boost::make_shared<Footprint>(getPeaks().getSchema());
993  int const xStart = getBBox().getMinX(), yStart = getBBox().getMinY();
994  std::vector<bool> rowBefore(width, false); // Representation of the previous row
995  std::vector<bool> rowNow(width, false); // Representation of this row
996  std::vector<bool> rowAfter(width, false); // Representation of the next row
997 
998  int yLast = yStart; // y value of last span we looked at
999  int const yEnd = _spans.back()->getY(); // y value of end span
1000 
1001  // Set rowNow, rowAfter
1002  SpanList::const_iterator readAhead = _spans.begin(); // Iterator for loading next row
1003  for (; readAhead != _spans.end() && (*readAhead)->getY() == yStart; ++readAhead) {
1004  std::fill(rowNow.begin() + (*readAhead)->getX0() - xStart,
1005  rowNow.begin() + (*readAhead)->getX1() + 1 - xStart,
1006  true);
1007  }
1008  for (; readAhead != _spans.end() && (*readAhead)->getY() == yStart + 1; ++readAhead) {
1009  std::fill(rowAfter.begin() + (*readAhead)->getX0() - xStart,
1010  rowAfter.begin() + (*readAhead)->getX1() + 1 - xStart,
1011  true);
1012  }
1013 
1014  for (SpanList::const_iterator ss = _spans.begin(); ss != _spans.end(); ++ss) {
1015  int const y = (*ss)->getY();
1016  if (y == yStart || y == yEnd) {
1017  // The whole span is on an edge
1018  edges->addSpanInSeries(y, (*ss)->getX0(), (*ss)->getX1());
1019  continue;
1020  }
1021  if (y != yLast) {
1022  // Move rows down
1023  rowBefore.assign(rowNow.begin(), rowNow.end());
1024  rowNow.assign(rowAfter.begin(), rowAfter.end());
1025  // Prepare the next row
1026  std::fill(rowAfter.begin(), rowAfter.end(), false);
1027  for (; readAhead != _spans.end() && (*readAhead)->getY() <= y; ++readAhead) {} // Moving only
1028  for (; readAhead != _spans.end() && (*readAhead)->getY() == y + 1; ++readAhead) {
1029  std::fill(rowAfter.begin() + (*readAhead)->getX0() - xStart,
1030  rowAfter.begin() + (*readAhead)->getX1() + 1 - xStart,
1031  true);
1032  }
1033  yLast = y;
1034  }
1035 
1036  // Look for edge in the current row
1037  int x0 = (*ss)->getX0();
1038  bool onEdge = true; // Are we on an edge? The first pixel is an edge
1039  for (int x = x0 + 1, i = x0 + 1 - xStart; x < (*ss)->getX1(); ++x, ++i) {
1040  if (onEdge) {
1041  if (rowBefore[i] && rowAfter[i]) {
1042  // We've come to the end of the edge
1043  onEdge = false;
1044  edges->addSpanInSeries(y, x0, x - 1);
1045  }
1046  } else if (!rowBefore[i] || !rowAfter[i]) {
1047  // We're on an edge now
1048  onEdge = true;
1049  x0 = x;
1050  }
1051  }
1052  // Last pixel is an edge
1053  int const x1 = (*ss)->getX1();
1054  if (onEdge) {
1055  edges->addSpanInSeries(y, x0, x1);
1056  } else {
1057  edges->addSpanInSeries(y, x1, x1);
1058  }
1059  }
1060  edges->normalize(); // Should be a no-op, but just in case...
1061 
1062  return edges;
1063 }
1064 
1065 
1071 bool _checkNormalized(Footprint const& foot) {
1072  Footprint copy(foot);
1073  copy.normalize();
1074  if (copy.getArea() != foot.getArea()) {
1075  return false;
1076  }
1077  if (copy.getSpans().size() != foot.getSpans().size()) {
1078  return false;
1079  }
1080  Footprint::SpanList const& spansa = foot.getSpans();
1081  Footprint::SpanList const& spansb = copy.getSpans();
1082  Footprint::SpanList::const_iterator spa = spansa.begin();
1083  Footprint::SpanList::const_iterator spb = spansb.begin();
1084  for (; spa != spansa.end(); ++spa, ++spb) {
1085  if ((*spa)->getY() != (*spb)->getY()) {
1086  return false;
1087  }
1088  if ((*spa)->getX0() != (*spb)->getX0()) {
1089  return false;
1090  }
1091  if ((*spa)->getX1() != (*spb)->getX1()) {
1092  return false;
1093  }
1094  }
1095  return true;
1096 }
1097 
1098 /************************************************************************************************************/
1099 
1100 template<typename MaskT>
1102  PTR(Footprint) const& fp,
1103  typename lsst::afw::image::Mask<MaskT>::Ptr const& mask,
1104  MaskT const bitmask
1105 ) {
1106  PTR(Footprint) newFp(new Footprint(fp->getPeaks().getSchema()));
1107  return newFp;
1108 }
1109 
1110 /************************************************************************************************************/
1111 
1112 template<typename MaskT>
1113 MaskT setMaskFromFootprint(
1114  image::Mask<MaskT> *mask,
1115  Footprint const& foot,
1116  MaskT const bitmask
1117 ) {
1118 
1119  int const width = static_cast<int>(mask->getWidth());
1120  int const height = static_cast<int>(mask->getHeight());
1121 
1122  for (Footprint::SpanList::const_iterator siter = foot.getSpans().begin();
1123  siter != foot.getSpans().end(); ++siter) {
1124  CONST_PTR(Span) span = *siter;
1125  int const y = span->getY() - mask->getY0();
1126  if (y < 0 || y >= height) {
1127  continue;
1128  }
1129 
1130  int x0 = span->getX0() - mask->getX0();
1131  int x1 = span->getX1() - mask->getX0();
1132  x0 = (x0 < 0) ? 0 : (x0 >= width ? width - 1 : x0);
1133  x1 = (x1 < 0) ? 0 : (x1 >= width ? width - 1 : x1);
1134 
1135  for (typename image::Image<MaskT>::x_iterator ptr = mask->x_at(x0, y),
1136  end = mask->x_at(x1 + 1, y); ptr != end; ++ptr) {
1137  *ptr |= bitmask;
1138  }
1139  }
1140 
1141  return bitmask;
1142 }
1143 
1144 /************************************************************************************************************/
1145 
1146 template<typename MaskT>
1148  image::Mask<MaskT> *mask,
1149  Footprint const& foot,
1150  MaskT const bitmask
1151 ) {
1152  int const width = static_cast<int>(mask->getWidth());
1153  int const height = static_cast<int>(mask->getHeight());
1154 
1155  for (Footprint::SpanList::const_iterator siter = foot.getSpans().begin();
1156  siter != foot.getSpans().end(); ++siter) {
1157  CONST_PTR(Span) span = *siter;
1158  int const y = span->getY() - mask->getY0();
1159  if (y < 0 || y >= height) {
1160  continue;
1161  }
1162 
1163  int x0 = span->getX0() - mask->getX0();
1164  int x1 = span->getX1() - mask->getX0();
1165  x0 = (x0 < 0) ? 0 : (x0 >= width ? width - 1 : x0);
1166  x1 = (x1 < 0) ? 0 : (x1 >= width ? width - 1 : x1);
1167 
1168  for (typename image::Image<MaskT>::x_iterator ptr = mask->x_at(x0, y),
1169  end = mask->x_at(x1 + 1, y); ptr != end; ++ptr) {
1170  *ptr &= ~bitmask;
1171  }
1172  }
1173 
1174  return bitmask;
1175 }
1176 
1177 /************************************************************************************************************/
1178 
1179 template<typename MaskT>
1181  image::Mask<MaskT> *mask,
1182  std::vector<PTR(Footprint)> const& footprints,
1183  MaskT const bitmask
1184 ) {
1185  for (std::vector<PTR(Footprint)>::const_iterator fiter = footprints.begin();
1186  fiter != footprints.end(); ++fiter) {
1187  (void)setMaskFromFootprint(mask, **fiter, bitmask);
1188  }
1189 
1190  return bitmask;
1191 }
1192 
1193 /************************************************************************************************************/
1194 
1195 template<typename MaskT>
1197  image::Mask<MaskT> *mask,
1198  CONST_PTR(std::vector<PTR(Footprint)>) const & footprints,
1199  MaskT const bitmask
1200  ) {
1201  return setMaskFromFootprintList(mask, *footprints, bitmask);
1202 }
1203 
1204 /************************************************************************************************************/
1205 namespace {
1206 template<typename ImageT>
1207 class SetFootprint : public FootprintFunctor<ImageT> {
1208 public:
1209  SetFootprint(ImageT const& image,
1210  typename ImageT::Pixel value) :
1211  FootprintFunctor<ImageT>(image), _value(value) {}
1212 
1213 
1214  void operator()(typename ImageT::xy_locator loc, int, int) {
1215  *loc = _value;
1216  }
1217 private:
1218  typename ImageT::Pixel _value;
1219 };
1220 }
1221 
1222 template<typename ImageT>
1223 typename ImageT::Pixel setImageFromFootprint(
1224  ImageT *image,
1225  Footprint const& foot,
1226  typename ImageT::Pixel const value
1227 ) {
1228  SetFootprint<ImageT> setit(*image, value);
1229  setit.apply(foot);
1230 
1231  return value;
1232 }
1233 
1234 template<typename ImageT>
1235 typename ImageT::Pixel setImageFromFootprintList(
1236  ImageT *image,
1237  CONST_PTR(std::vector<PTR(Footprint)>) footprints,
1238  typename ImageT::Pixel const value
1239  ) {
1240  return setImageFromFootprintList(image, *footprints, value);
1241 }
1242 
1243 template<typename ImageT>
1244 typename ImageT::Pixel setImageFromFootprintList(
1245  ImageT *image,
1246  std::vector<PTR(Footprint)> const& footprints,
1247  typename ImageT::Pixel const value
1248 ) {
1249  SetFootprint<ImageT> setit(*image, value);
1250  for (std::vector<PTR(Footprint)>::const_iterator fiter = footprints.begin(),
1251  end = footprints.end(); fiter != end; ++fiter) {
1252  setit.apply(**fiter);
1253  }
1254 
1255  return value;
1256 }
1257 
1258 /************************************************************************************************************/
1259 /*
1260  * Worker routine for the pmSetFootprintArrayIDs/pmSetFootprintID (and pmMergeFootprintArrays)
1261  */
1262 template <typename IDPixelT>
1263 static void set_footprint_id(
1264  typename image::Image<IDPixelT>::Ptr idImage, // the image to set
1265  Footprint const& foot, // the footprint to insert
1266  int const id, // the desired ID
1267  int dx=0, int dy=0 // Add these to all x/y in the Footprint
1268 ) {
1269  for (Footprint::SpanList::const_iterator i = foot.getSpans().begin();
1270  i != foot.getSpans().end(); ++i) {
1271  CONST_PTR(Span) span = *i;
1272  for (typename image::Image<IDPixelT>::x_iterator ptr =
1273  idImage->x_at(span->getX0() + dx, span->getY() + dy),
1274  end = ptr + span->getWidth(); ptr != end; ++ptr) {
1275  *ptr = id;
1276  }
1277  }
1278 }
1279 
1280 template <typename IDPixelT>
1281 static void
1282 set_footprint_array_ids(
1283  typename image::Image<IDPixelT>::Ptr idImage, // the image to set
1284  std::vector<PTR(Footprint)> const& footprints, // the footprints to insert
1285  bool const relativeIDs // show IDs starting at 0, not Footprint->id
1286 ) {
1287  int id = 0; // first index will be 1
1288 
1289  for (std::vector<PTR(Footprint)>::const_iterator fiter = footprints.begin();
1290  fiter != footprints.end(); ++fiter) {
1291  CONST_PTR(Footprint) foot = *fiter;
1292 
1293  if (relativeIDs) {
1294  ++id;
1295  } else {
1296  id = foot->getId();
1297  }
1298 
1299  set_footprint_id<IDPixelT>(idImage, *foot, id);
1300  }
1301 }
1302 
1303 template void set_footprint_array_ids<int>(
1304  image::Image<int>::Ptr idImage,
1305  std::vector<PTR(Footprint)> const& footprints,
1306  bool const relativeIDs);
1307 
1308 /******************************************************************************/
1309 /*
1310  * Set an image to the value of footprint's ID wherever they may fall
1311  *
1312  * @param footprints the footprints to insert
1313  * @param relativeIDs show the IDs starting at 1, not pmFootprint->id
1314  */
1315 template <typename IDImageT>
1316 typename boost::shared_ptr<image::Image<IDImageT> > setFootprintArrayIDs(
1317  std::vector<PTR(Footprint)> const& footprints,
1318  bool const relativeIDs
1319 ) {
1320  std::vector<PTR(Footprint)>::const_iterator fiter = footprints.begin();
1321  if (fiter == footprints.end()) {
1322  throw LSST_EXCEPT(
1323  lsst::pex::exceptions::InvalidParameterError,
1324  "You didn't provide any footprints"
1325  );
1326  }
1327  CONST_PTR(Footprint) foot = *fiter;
1328 
1329  typename image::Image<IDImageT>::Ptr idImage(
1330  new image::Image<IDImageT>(foot->getRegion())
1331  );
1332  *idImage = 0;
1333  /*
1334  * do the work
1335  */
1336  set_footprint_array_ids<IDImageT>(idImage, footprints, relativeIDs);
1337 
1338  return idImage;
1339 }
1340 
1341 template image::Image<int>::Ptr setFootprintArrayIDs(
1342  std::vector<PTR(Footprint)> const& footprints,
1343  bool const relativeIDs);
1344 /*
1345  * Set an image to the value of Footprint's ID wherever it may fall
1346  */
1347 template <typename IDImageT>
1348 typename boost::shared_ptr<image::Image<IDImageT> > setFootprintID(
1349  CONST_PTR(Footprint)& foot, // the Footprint to insert
1350  int const id // the desired ID
1351  ) {
1352  typename image::Image<IDImageT>::Ptr idImage(new image::Image<IDImageT>(foot->getBBox()));
1353  *idImage = 0;
1354  /*
1355  * do the work
1356  */
1357  set_footprint_id<IDImageT>(idImage, *foot, id);
1358 
1359  return idImage;
1360 }
1361 
1362 template image::Image<int>::Ptr setFootprintID(CONST_PTR(Footprint)& foot, int const id);
1363 
1364 namespace {
1371 class StructuringElement
1372 {
1373 public:
1374  enum class Shape { CIRCLE, DIAMOND };
1375  typedef std::vector<Span>::const_iterator const_iterator;
1376  StructuringElement(Shape shape, int radius);
1377  StructuringElement(int left, int right, int up, int down);
1378  const_iterator begin() const { return _widths.begin(); }
1379  const_iterator end() const { return _widths.end(); }
1380  int getYRange() const { return _yRange; }
1381 
1382 private:
1383  std::vector<Span> _widths;
1384  int _yRange;
1385 };
1386 
1392 StructuringElement::StructuringElement(Shape shape, int radius) {
1393  _yRange = 2*radius + 1;
1394  _widths.reserve(_yRange);
1395  switch (shape) {
1396  case Shape::CIRCLE:
1397  for (auto dy = -radius; dy <= radius; dy++) {
1398  int dx = static_cast<int>(sqrt(radius*radius - dy*dy));
1399  _widths.push_back(Span(dy, -dx, dx));
1400  }
1401  break;
1402  case Shape::DIAMOND:
1403  for (auto dy = -radius; dy <= radius; dy++) {
1404  int dx = radius - abs(dy);
1405  _widths.push_back(Span(dy, -dx, dx));
1406  }
1407  break;
1408  }
1409 }
1410 
1415 StructuringElement::StructuringElement(int left, int right, int up, int down) {
1416  _yRange = up + down + 1;
1417  _widths.reserve(_yRange);
1418  for (auto dy = 1; dy <= up; dy++) {
1419  _widths.push_back(Span(dy, 0, 0));
1420  }
1421  for (auto dy = -1; dy >= -down; dy--) {
1422  _widths.push_back(Span(dy, 0, 0));
1423  }
1424  _widths.push_back(Span(0, -left, right));
1425 }
1426 
1431 PTR(Footprint) growFootprintImpl(
1432  Footprint const& foot,
1433  StructuringElement const& element
1434 ) {
1435  // Create an empty footprint covering foot's region.
1436  PTR(Footprint) grown(new Footprint(0, foot.getRegion()));
1437 
1438  // We use a map of (y coordinate) to set of (xmin, xmax) pairs to describe
1439  // the spans being constructed. In this way, we ensure that the spans are
1440  // always sorted by increasing y, xmin.
1441  std::map<int, std::set<std::pair<int, int>>> spans;
1442 
1443  // Iterate over foot & structuring element building up a collection of
1444  // spans which should be added to the footprint.
1445  for (auto spanIter = foot.getSpans().cbegin(); spanIter != foot.getSpans().cend(); ++spanIter) {
1446  for (auto elementIter = element.begin(); elementIter != element.end(); ++elementIter) {
1447  int const xmin = (*spanIter)->getX0() + elementIter->getX0();
1448  int const xmax = (*spanIter)->getX1() + elementIter->getX1();
1449  int const yval = (*spanIter)->getY() + elementIter->getY();
1450  spans[yval].insert(std::make_pair(xmin, xmax));
1451 
1452  // Merge overlapping spans at this y coordinate, thereby ensuring
1453  // that the proto-footprint remains normalized.
1454  std::set<std::pair<int, int>> newSpans;
1455  for (auto span = spans[yval].cbegin(); span != spans[yval].cend(); ++span) {
1456  int const start = span->first;
1457  int end = span->second;
1458  // Check for end+1 because end value is inclusive. That is, if
1459  // one span terminates at x=N and another begins at x=N+1,
1460  // those spans are contiguous.
1461  while (span != spans[yval].cend() && span->first <= (end+1)) {
1462  end = std::max(end, span++->second);
1463  }
1464  newSpans.insert(std::make_pair(start, end));
1465  --span; // "rewind" to the last span included
1466  }
1467  std::swap(spans[yval], newSpans);
1468  }
1469  }
1470 
1471  // Now append the spans to the output footprint, making use of the fact
1472  // that they are already normalized.
1473  for (auto y = spans.cbegin(); y != spans.cend(); ++y) {
1474  for (auto x = y->second.cbegin(); x != y->second.cend(); ++x) {
1475  grown->addSpanInSeries(y->first, x->first, x->second);
1476  }
1477  }
1478 
1479  // Copy over peaks from the original footprint
1480  grown->getPeaks() = PeakCatalog(foot.getPeaks().getTable(), foot.getPeaks().begin(),
1481  foot.getPeaks().end(), true);
1482 
1483  return grown;
1484 }
1485 
1495 struct PrimaryRun {
1496  int m, y, xmin, xmax;
1497 };
1498 
1503 bool comparePrimaryRun(PrimaryRun const& first, PrimaryRun const& second) {
1504  if (first.y != second.y) {
1505  return first.y < second.y;
1506  } else if (first.m != second.m) {
1507  return first.m < second.m;
1508  } else {
1509  return first.xmin < second.xmin;
1510  }
1511 }
1512 
1513 class ComparePrimaryRunY{
1514 public:
1515  bool operator()(PrimaryRun const& pr, int yval) {
1516  return pr.y < yval;
1517  }
1518  bool operator()(int yval, PrimaryRun const& pr) {
1519  return yval < pr.y;
1520  }
1521 };
1522 
1523 class ComparePrimaryRunM{
1524 public:
1525  bool operator()(PrimaryRun const& pr, int mval) {
1526  return pr.m < mval;
1527  }
1528  bool operator()(int mval, PrimaryRun const& pr) {
1529  return mval < pr.m;
1530  }
1531 };
1532 
1537 PTR(Footprint) shrinkFootprintImpl(
1538  Footprint const& foot,
1539  StructuringElement const& element
1540 ) {
1541  // Create an empty FootprintSet covering the input region
1542  PTR(Footprint) shrunk(new Footprint(0, foot.getRegion()));
1543 
1544  // Calculate all possible primary runs.
1545  std::vector<PrimaryRun> primaryRuns;
1546  for (auto spanIter = foot.getSpans().begin(); spanIter != foot.getSpans().end(); ++spanIter) {
1547  int m = 0;
1548  for (auto it = element.begin(); it != element.end(); ++it, ++m) {
1549  if ((it->getX1() - it->getX0()) <= ((*spanIter)->getX1() - (*spanIter)->getX0())) {
1550  int xmin = (*spanIter)->getX0() - it->getX0();
1551  int xmax = (*spanIter)->getX1() - it->getX1();
1552  int y = (*spanIter)->getY() - it->getY();
1553  primaryRuns.push_back(PrimaryRun({m, y, xmin, xmax}));
1554  }
1555  }
1556  }
1557 
1558  // Iterate over the primary runs in such a way that we consider all values
1559  // of m for given y, then all m for y+1, etc.
1560  std::sort(primaryRuns.begin(), primaryRuns.end(), comparePrimaryRun);
1561 
1562  for (int y = primaryRuns.front().y; y <= primaryRuns.back().y; ++y) {
1563  auto yRange = std::equal_range(primaryRuns.begin(), primaryRuns.end(), y, ComparePrimaryRunY());
1564 
1565  // Discard runs for any value of y for which we find fewer groups than
1566  // M, the total Y range of the structuring element. This is step 3.1
1567  // of the Kim et al. algorithm.
1568  if (std::distance(yRange.first, yRange.second) < element.getYRange()) {
1569  continue;
1570  }
1571 
1572  // "good" runs are those which are covered by each value of m, ie by
1573  // each row in the structuring element. Our algorithm will consider
1574  // each value of m in turn, gradually whittling down the list of good
1575  // runs, then finally convert the remainder into Spans and add them to
1576  // the shrunken Footprint.
1577  std::list<PrimaryRun> goodRuns;
1578 
1579  for (int m = 0; m < element.getYRange(); ++m) {
1580  auto mRange = std::equal_range(yRange.first, yRange.second, m, ComparePrimaryRunM());
1581  if (mRange.first == mRange.second) {
1582  // If a particular m is missing, we know that this y contains
1583  // no good runs; this is equivalent to Kim et al. step 3.2.
1584  goodRuns.clear();
1585  } else {
1586  // Consolidate all primary runs at this m so that they
1587  // don't overlap.
1588  std::list<PrimaryRun> candidateRuns;
1589  int start_x = mRange.first->xmin;
1590  int end_x = mRange.first->xmax;
1591  for (auto run = mRange.first+1; run != mRange.second; ++run) {
1592  if (run->xmin > end_x) {
1593  // Start of a new run
1594  candidateRuns.push_back(PrimaryRun{m, y, start_x, end_x});
1595  start_x = run->xmin;
1596  end_x = run->xmax;
1597  } else {
1598  // Continuation of an existing run
1599  end_x = run->xmax;
1600  }
1601  }
1602  candidateRuns.push_back(PrimaryRun{m, y, start_x, end_x});
1603 
1604  // Otherwise, calculate the intersection of candidate runs at
1605  // this m with good runs from all previous m.
1606  if (m == 0) {
1607  // For m = 0 we have nothing to compare to; all runs are
1608  // accepted.
1609  std::swap(goodRuns, candidateRuns);
1610  } else {
1611  std::list<PrimaryRun> newlist;
1612  for (auto good = goodRuns.begin(); good != goodRuns.end(); ++good) {
1613  for (auto cand = candidateRuns.begin(); cand != candidateRuns.end(); ++cand) {
1614  int start = std::max(good->xmin, cand->xmin);
1615  int end = std::min(good->xmax, cand->xmax);
1616  if (end >= start) {
1617  newlist.push_back(PrimaryRun({m, y, start, end}));
1618  }
1619  }
1620  }
1621  std::swap(newlist, goodRuns);
1622  }
1623  }
1624  }
1625  for (auto run = goodRuns.begin(); run != goodRuns.end(); ++run) {
1626  shrunk->addSpan(run->y, run->xmin, run->xmax);
1627  }
1628  }
1629 
1630  shrunk->normalize();
1631 
1632  // Peaks from the original footprint have not yet been added to the shrunken footprint.
1633  // Iterate over peaks from the original footprint and add them IF they are contained
1634  // within the shrunken footprint.
1635  for (auto peakIter = foot.getPeaks().begin(); peakIter != foot.getPeaks().end(); peakIter++) {
1636  if (shrunk->contains(peakIter->getI())) {
1637  shrunk->getPeaks().addNew()->assign(*peakIter);
1638  }
1639  }
1640  return shrunk;
1641 }
1642 }
1643 
1644 /************************************************************************************************************/
1645 
1646 namespace {
1647  PTR(Footprint) _mergeFootprints(Footprint const& aFoot, Footprint const& bFoot) {
1648  PTR(Footprint) foot(new Footprint());
1649 
1650  const PeakCatalog& aPeak = aFoot.getPeaks();
1651  const PeakCatalog& bPeak = bFoot.getPeaks();
1652  PeakCatalog& peaks = foot->getPeaks();
1653  if (aPeak.empty()) {
1654  if (!bPeak.empty()) {
1655  peaks = PeakCatalog(bPeak.getTable(), bPeak.begin(), bPeak.end(), true);
1656  }
1657  } else {
1658  if (bPeak.empty()) {
1659  peaks = PeakCatalog(aPeak.getTable(), aPeak.begin(), aPeak.end(), true);
1660  } else {
1661  if (aPeak.getSchema() == bPeak.getSchema()) {
1662  // use schema A, as it's the same as schema B
1663  peaks = PeakCatalog(aPeak.getTable());
1664  peaks.reserve(aPeak.size() + bPeak.size());
1665  peaks.insert(peaks.end(), aPeak.begin(), aPeak.end(), true);
1666  peaks.insert(peaks.end(), bPeak.begin(), bPeak.end(), true);
1667  } else {
1668  throw LSST_EXCEPT(
1669  pex::exceptions::InvalidParameterError,
1670  "Cannot merge Footprints when Peaks have different Schemas"
1671  );
1672  }
1673  }
1674  }
1675 
1676  Footprint::SpanList const& aSpans = aFoot.getSpans();
1677  Footprint::SpanList const& bSpans = bFoot.getSpans();
1678  Footprint::SpanList::const_iterator aSpan = aSpans.begin();
1679  Footprint::SpanList::const_iterator bSpan = bSpans.begin();
1680  Footprint::SpanList::const_iterator aEnd = aSpans.end();
1681  Footprint::SpanList::const_iterator bEnd = bSpans.end();
1682 
1683  foot->getSpans().reserve(std::max(aSpans.size(), bSpans.size()));
1684 
1685  while ((aSpan != aEnd) && (bSpan != bEnd)) {
1686  int y = (*aSpan)->getY();
1687  int x0 = (*aSpan)->getX0();
1688  int x1 = (*aSpan)->getX1();
1689  int yb = (*bSpan)->getY();
1690  int xb0 = (*bSpan)->getX0();
1691  int xb1 = (*bSpan)->getX1();
1692 
1693  if ((y < yb) || (y == yb && (x1 < (xb0-1)))) {
1694  // A is earlier -- add A
1695  foot->addSpanInSeries(y, x0, x1);
1696  ++aSpan;
1697  continue;
1698  }
1699  if ((yb < y) || (y == yb && (xb1 < (x0-1)))) {
1700  // B is earlier -- add B
1701  foot->addSpanInSeries(yb, xb0, xb1);
1702  ++bSpan;
1703  continue;
1704  }
1705 
1706  assert(yb == y);
1707  // Overlap -- find connected spans from both iterators.
1708  x0 = std::min(x0, xb0);
1709  x1 = std::max(x1, xb1);
1710  // Union all connected spans
1711  ++aSpan;
1712  ++bSpan;
1713  while (true) {
1714  if ((aSpan != aEnd) &&
1715  ((*aSpan)->getY() == y) &&
1716  ((*aSpan)->getX0() <= (x1+1))) {
1717  // *aSpan continues this span.
1718  x1 = std::max(x1, (*aSpan)->getX1());
1719  ++aSpan;
1720  continue;
1721  }
1722  if ((bSpan != bEnd) &&
1723  ((*bSpan)->getY() == y) &&
1724  ((*bSpan)->getX0() <= (x1+1))) {
1725  // *bSpan continues this span.
1726  x1 = std::max(x1, (*bSpan)->getX1());
1727  ++bSpan;
1728  continue;
1729  }
1730  break;
1731  }
1732  foot->addSpanInSeries(y, x0, x1);
1733  }
1734  // At this point either "aSpan" or "bSpan" is at the end.
1735 
1736  // Add any remaining spans from "A".
1737  for (; aSpan != aEnd; ++aSpan) {
1738  foot->addSpanInSeries((*aSpan)->getY(), (*aSpan)->getX0(), (*aSpan)->getX1());
1739  }
1740  // Add any remaining spans from "B".
1741  for (; bSpan != bEnd; ++bSpan) {
1742  foot->addSpanInSeries((*bSpan)->getY(), (*bSpan)->getX0(), (*bSpan)->getX1());
1743  }
1744  return foot;
1745  }
1746 }
1747 
1748 PTR(Footprint) mergeFootprints(Footprint& foot1, Footprint& foot2) {
1749  foot1.normalize();
1750  foot2.normalize();
1751  return _mergeFootprints(foot1, foot2);
1752 }
1753 
1754 PTR(Footprint) mergeFootprints(Footprint const& foot1, Footprint const& foot2) {
1755  if (!foot1.isNormalized() || !foot2.isNormalized()) {
1756  throw LSST_EXCEPT(
1757  lsst::pex::exceptions::InvalidParameterError,
1758  "mergeFootprints(const Footprints) requires normalize()d Footprints.");
1759  }
1760  return _mergeFootprints(foot1, foot2);
1761 }
1762 
1763 /************************************************************************************************************/
1764 
1765 void nearestFootprint(std::vector<PTR(Footprint)> const& foots,
1768 {
1769  /*
1770  * insert the footprints into an image, set all the pixels to the
1771  * Manhattan distance from the nearest set pixel.
1772  */
1773  typedef boost::uint16_t dtype;
1774  typedef boost::uint16_t itype;
1775 
1776  const itype nil = 0xffff;
1777 
1778  const geom::Box2I bbox = argmin->getBBox();
1779  *argmin = 0;
1780  *dist = 0;
1781 
1782  int const x0 = bbox.getMinX();
1783  int const y0 = bbox.getMinY();
1784 
1785  for (size_t i=0; i<foots.size(); ++i) {
1786  // Set all the pixels in the footprint to 1
1787  set_footprint_id<itype>(argmin, *foots[i], i, -x0, -y0);
1788  set_footprint_id<dtype>(dist, *foots[i], 1, -x0, -y0);
1789  }
1790 
1791  int const height = dist->getHeight();
1792  int const width = dist->getWidth();
1793 
1794  // traverse from bottom left to top right
1795  for (int y = 0; y != height; ++y) {
1796  image::Image<dtype>::xy_locator dim = dist->xy_at(0, y);
1797  image::Image<itype>::xy_locator aim = argmin->xy_at(0, y);
1798  for (int x = 0; x != width; ++x, ++dim.x(), ++aim.x()) {
1799  if (dim(0, 0) == 1) {
1800  // first pass and pixel was on, it gets a zero
1801  dim(0, 0) = 0;
1802  // its argmin is already set
1803  } else {
1804  // pixel was off. It is at most the sum of lengths of
1805  // the array away from a pixel that is on
1806  dim(0, 0) = width + height;
1807  aim(0, 0) = nil;
1808  // or one more than the pixel to the north
1809  if (y > 0) {
1810  dtype ndist = dim(0,-1) + 1;
1811  if (ndist < dim(0,0)) {
1812  dim(0,0) = ndist;
1813  aim(0,0) = aim(0,-1);
1814  }
1815  }
1816  // or one more than the pixel to the west
1817  if (x > 0) {
1818  dtype ndist = dim(-1,0) + 1;
1819  if (ndist < dim(0,0)) {
1820  dim(0,0) = ndist;
1821  aim(0,0) = aim(-1,0);
1822  }
1823  }
1824  }
1825  }
1826  }
1827  // traverse from top right to bottom left
1828  for (int y = height - 1; y >= 0; --y) {
1829  image::Image<dtype>::xy_locator dim = dist->xy_at(width-1, y);
1830  image::Image<itype>::xy_locator aim = argmin->xy_at(width-1, y);
1831  for (int x = width - 1; x >= 0; --x, --dim.x(), --aim.x()) {
1832  // either what we had on the first pass or one more than
1833  // the pixel to the south
1834  if (y + 1 < height) {
1835  dtype ndist = dim(0,1) + 1;
1836  if (ndist < dim(0,0)) {
1837  dim(0,0) = ndist;
1838  aim(0,0) = aim(0,1);
1839  }
1840  }
1841  // or one more than the pixel to the east
1842  if (x + 1 < width) {
1843  dtype ndist = dim(1,0) + 1;
1844  if (ndist < dim(0,0)) {
1845  dim(0,0) = ndist;
1846  aim(0,0) = aim(1,0);
1847  }
1848  }
1849  }
1850  }
1851 }
1852 
1853 PTR(Footprint) growFootprint(
1854  Footprint const& foot,
1855  int nGrow,
1856  bool isotropic
1857 ) {
1858  if (nGrow <= 0 || foot.getNpix() == 0 ) {
1859  // Return a new footprint equal to the input.
1860  return PTR(Footprint)(new Footprint(foot));
1861  }
1862 
1863  // An isotropic grow is equivalent to growing with a circular structuring
1864  // element, while a Manhattan grow is equivalent to growing with a
1865  // diamond-shaped element.
1867  Shape shape = isotropic ? Shape::CIRCLE : Shape::DIAMOND;
1868  return growFootprintImpl(foot, StructuringElement(shape, nGrow));
1869 }
1870 
1871 PTR(Footprint) growFootprint(PTR(Footprint) const& foot, int nGrow, bool isotropic) {
1872  return growFootprint(*foot, nGrow, isotropic);
1873 }
1874 
1875 PTR(Footprint) growFootprint(Footprint const& foot,
1876  int nGrow,
1877  bool left,
1878  bool right,
1879  bool up,
1880  bool down
1881  )
1882 {
1883  if (nGrow <= 0 || foot.getNpix() == 0 ) {
1884  // Return a new footprint equal to the input.
1885  return PTR(Footprint)(new Footprint(foot));
1886  }
1887  return growFootprintImpl(foot, StructuringElement(left ? nGrow: 0, right ? nGrow : 0,
1888  up ? nGrow : 0, down ? nGrow : 0));
1889 }
1890 
1891 
1892 PTR(Footprint) shrinkFootprint(
1893  Footprint const& foot,
1894  int nShrink,
1895  bool isotropic
1896 ) {
1898  Shape shape = isotropic ? Shape::CIRCLE : Shape::DIAMOND;
1899  return shrinkFootprintImpl(foot, StructuringElement(shape, nShrink));
1900 }
1901 
1902 /************************************************************************************************************/
1903 
1904 std::vector<geom::Box2I> footprintToBBoxList(Footprint const& foot) {
1905  typedef boost::uint16_t ImageT;
1906  geom::Box2I fpBBox = foot.getBBox();
1907  image::Image<ImageT>::Ptr idImage(
1908  new image::Image<ImageT>(fpBBox.getDimensions())
1909  );
1910  *idImage = 0;
1911  int const height = fpBBox.getHeight();
1912  geom::Extent2I shift(fpBBox.getMinX(), fpBBox.getMinY());
1913  foot.insertIntoImage(*idImage, 1, fpBBox);
1914 
1915  std::vector<geom::Box2I> bboxes;
1916  /*
1917  * Our strategy is to find a row of pixels in the Footprint and interpret it as the first
1918  * row of a rectangular set of pixels. We then extend this rectangle upwards as far as it
1919  * will go, and define that as a BBox. We clear all those pixels, and repeat until there
1920  * are none left. I.e. a Footprint will get cut up like this:
1921  *
1922  * .555...
1923  * 22.3314
1924  * 22.331.
1925  * .000.1.
1926  * (as shown in Footprint_1.py)
1927  */
1928 
1929  int y0 = 0; // the first row with non-zero pixels in it
1930  while (y0 < height) {
1931  geom::Box2I bbox; // our next BBox
1932  for (int y = y0; y != height; ++y) {
1933  // Look for a set pixel in this row
1934  image::Image<ImageT>::x_iterator begin = idImage->row_begin(y), end = idImage->row_end(y);
1935  image::Image<ImageT>::x_iterator first = std::find(begin, end, 1);
1936 
1937  if (first != end) { // A pixel is set in this row
1938  image::Image<ImageT>::x_iterator last = std::find(first, end, 0) - 1;
1939  int const x0 = first - begin;
1940  int const x1 = last - begin;
1941 
1942  std::fill(first, last + 1, 0); // clear pixels; we don't want to see them again
1943 
1944  bbox.include(geom::Point2I(x0, y)); // the LLC
1945  bbox.include(geom::Point2I(x1, y)); // the LRC; initial guess for URC
1946 
1947  // we found at least one pixel so extend the BBox upwards
1948  for (++y; y != height; ++y) {
1949  if (std::find(idImage->at(x0, y), idImage->at(x1 + 1, y), 0) != idImage->at(x1 + 1, y)) {
1950  break; // some pixels weren't set, so the BBox stops here, (actually in previous row)
1951  }
1952  std::fill(idImage->at(x0, y), idImage->at(x1 + 1, y), 0);
1953 
1954  bbox.include(geom::Point2I(x1, y)); // the new URC
1955  }
1956 
1957  bbox.shift(shift);
1958  bboxes.push_back(bbox);
1959  } else {
1960  y0 = y + 1;
1961  }
1962  break;
1963  }
1964  }
1965 
1966  return bboxes;
1967 }
1968 
1969 
1970 template<typename ImageOrMaskedImageT>
1971 void
1972 copyWithinFootprint(Footprint const& foot,
1973  PTR(ImageOrMaskedImageT) const input,
1974  PTR(ImageOrMaskedImageT) output) {
1975  Footprint::SpanList spans = foot.getSpans();
1976 
1977  int const inX0 = input->getX0(), inY0 = input->getY0();
1978  int const outX0 = output->getX0(), outY0 = output->getY0();
1979 
1980  int const xMin = std::max(inX0, outX0);
1981  int const xMax = std::min(input->getWidth() + inX0, output->getWidth() + outX0) - 1;
1982  for (Footprint::SpanList::iterator sp = spans.begin();
1983  sp != spans.end(); ++sp) {
1984  int const y = (*sp)->getY();
1985  int const x0 = (*sp)->getX0();
1986  int const x1 = (*sp)->getX1();
1987 
1988  int const yInput = y - inY0, yOutput = y - outY0;
1989  if (yInput < 0 || yInput >= input->getHeight() || yOutput < 0 || yOutput >= output->getHeight()) {
1990  continue;
1991  }
1992 
1993  int const xStart = std::max(x0, xMin); // Starting position in x, parent frame
1994  int const xStop = std::min(x1, xMax); // Stopping position (inclusive) in x, parent frame
1995 
1996  typename ImageOrMaskedImageT::const_x_iterator initer = input->x_at(xStart - inX0, yInput);
1997  typename ImageOrMaskedImageT::x_iterator outiter = output->x_at(xStart - outX0, yOutput);
1998  for (int x = xStart; x <= xStop; ++x, ++initer, ++outiter) {
1999  *outiter = *initer;
2000  }
2001  }
2002 }
2003 
2004 
2005 
2006 
2007 #if 0
2008 
2009 /************************************************************************************************************/
2010 /*
2011  * Grow a psArray of pmFootprints isotropically by r pixels, returning a new psArray of new pmFootprints
2012  */
2013 psArray *pmGrowFootprintArray(psArray const *footprints, // footprints to grow
2014  int r) { // how much to grow each footprint
2015  assert (footprints->n == 0 || pmIsFootprint(footprints->data[0]));
2016 
2017  if (footprints->n == 0) { // we don't know the size of the footprint's region
2018  return psArrayAlloc(0);
2019  }
2020  /*
2021  * We'll insert the footprints into an image, then convolve with a disk,
2022  * then extract a footprint from the result --- this is magically what we want.
2023  */
2024  psImage *idImage = pmSetFootprintArrayIDs(footprints, true);
2025  if (r <= 0) {
2026  r = 1; // r == 1 => no grow
2027  }
2028  psKernel *circle = psKernelAlloc(-r, r, -r, r);
2029  assert (circle->image->numRows == 2*r + 1 && circle->image->numCols == circle->image->numRows);
2030  for (int i = 0; i <= r; ++i) {
2031  for (int j = 0; j <= r; ++j) {
2032  if (i*i + j*j <= r*r) {
2033  circle->kernel[i][j] =
2034  circle->kernel[i][-j] =
2035  circle->kernel[-i][j] =
2036  circle->kernel[-i][-j] = 1;
2037  }
2038  }
2039  }
2040 
2041  psImage *grownIdImage = psImageConvolveDirect(idImage, circle); // Here's the actual grow step
2042  psFree(circle);
2043 
2044  psArray *grown = pmFindFootprints(grownIdImage, 0.5, 1); // and here we rebuild the grown footprints
2045  assert (grown != NULL);
2046  psFree(idImage);
2047  psFree(grownIdImage);
2048  /*
2049  * Now assign the peaks appropriately. We could do this more efficiently
2050  * using grownIdImage (which we just freed), but this is easy and probably fast enough
2051  */
2052  psArray const *peaks = pmFootprintArrayToPeaks(footprints);
2053  pmPeaksAssignToFootprints(grown, peaks);
2054  psFree((psArray *)peaks);
2055 
2056  return grown;
2057 }
2058 
2059 /************************************************************************************************************/
2060 /*
2061  * Merge together two psArrays of pmFootprints neither of which is damaged.
2062  *
2063  * The returned psArray may contain elements of the inital psArrays (with
2064  * their reference counters suitable incremented)
2065  */
2066 psArray *pmMergeFootprintArrays(psArray const *footprints1, // one set of footprints
2067  psArray const *footprints2, // the other set
2068  int const includePeaks) { // which peaks to set? 0x1 => footprints1, 0x2 => 2
2069  assert (footprints1->n == 0 || pmIsFootprint(footprints1->data[0]));
2070  assert (footprints2->n == 0 || pmIsFootprint(footprints2->data[0]));
2071 
2072  if (footprints1->n == 0 || footprints2->n == 0) { // nothing to do but put copies on merged
2073  psArray const *old = (footprints1->n == 0) ? footprints2 : footprints1;
2074 
2075  psArray *merged = psArrayAllocEmpty(old->n);
2076  for (int i = 0; i < old->n; ++i) {
2077  psArrayAdd(merged, 1, old->data[i]);
2078  }
2079 
2080  return merged;
2081  }
2082  /*
2083  * We have real work to do as some pmFootprints in footprints2 may overlap
2084  * with footprints1
2085  */
2086  {
2087  pmFootprint *fp1 = footprints1->data[0];
2088  pmFootprint *fp2 = footprints2->data[0];
2089  if (fp1->region.x0 != fp2->region.x0 ||
2090  fp1->region.x1 != fp2->region.x1 ||
2091  fp1->region.y0 != fp2->region.y0 ||
2092  fp1->region.y1 != fp2->region.y1) {
2093  psError(PS_ERR_BAD_PARAMETER_SIZE, true,
2094  "The two pmFootprint arrays correspnond to different-sized regions");
2095  return NULL;
2096  }
2097  }
2098  /*
2099  * We'll insert first one set of footprints then the other into an image, then
2100  * extract a footprint from the result --- this is magically what we want.
2101  */
2102  psImage *idImage = pmSetFootprintArrayIDs(footprints1, true);
2103  set_footprint_array_ids(idImage, footprints2, true);
2104 
2105  psArray *merged = pmFindFootprints(idImage, 0.5, 1);
2106  assert (merged != NULL);
2107  psFree(idImage);
2108  /*
2109  * Now assign the peaks appropriately. We could do this more efficiently
2110  * using idImage (which we just freed), but this is easy and probably fast enough
2111  */
2112  if (includePeaks & 0x1) {
2113  psArray const *peaks = pmFootprintArrayToPeaks(footprints1);
2114  pmPeaksAssignToFootprints(merged, peaks);
2115  psFree((psArray *)peaks);
2116  }
2117 
2118  if (includePeaks & 0x2) {
2119  psArray const *peaks = pmFootprintArrayToPeaks(footprints2);
2120  pmPeaksAssignToFootprints(merged, peaks);
2121  psFree((psArray *)peaks);
2122  }
2123 
2124  return merged;
2125 }
2126 
2127 /************************************************************************************************************/
2128 /*
2129  * Given a psArray of pmFootprints and another of pmPeaks, assign the peaks to the
2130  * footprints in which that fall; if they _don't_ fall in a footprint, add a suitable
2131  * one to the list.
2132  */
2133 psErrorCode
2134 pmPeaksAssignToFootprints(psArray *footprints, // the pmFootprints
2135  psArray const *peaks) { // the pmPeaks
2136  assert (footprints != NULL);
2137  assert (footprints->n == 0 || pmIsFootprint(footprints->data[0]));
2138  assert (peaks != NULL);
2139  assert (peaks->n == 0 || pmIsPeak(peaks->data[0]));
2140 
2141  if (footprints->n == 0) {
2142  if (peaks->n > 0) {
2143  return psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Your list of footprints is empty");
2144  }
2145  return PS_ERR_NONE;
2146  }
2147  /*
2148  * Create an image filled with the object IDs, and use it to assign pmPeaks to the
2149  * objects
2150  */
2151  psImage *ids = pmSetFootprintArrayIDs(footprints, true);
2152  assert (ids != NULL);
2153  assert (ids->type.type == PS_TYPE_S32);
2154  int const y0 = ids->y0;
2155  int const x0 = ids->x0;
2156  int const numRows = ids->numRows;
2157  int const numCols = ids->numCols;
2158 
2159  for (int i = 0; i < peaks->n; ++i) {
2160  pmPeak *peak = peaks->data[i];
2161  int const x = peak->x - x0;
2162  int const y = peak->y - y0;
2163 
2164  assert (x >= 0 && x < numCols && y >= 0 && y < numRows);
2165  int id = ids->data.S32[y][x - x0];
2166 
2167  if (id == 0) { // peak isn't in a footprint, so make one for it
2168  pmFootprint *nfp = pmFootprintAlloc(1, ids);
2169  pmFootprintAddSpan(nfp, y, x, x);
2170  psArrayAdd(footprints, 1, nfp);
2171  psFree(nfp);
2172  id = footprints->n;
2173  }
2174 
2175  assert (id >= 1 && id <= footprints->n);
2176  pmFootprint *fp = footprints->data[id - 1];
2177  psArrayAdd(fp->peaks, 5, peak);
2178  }
2179 
2180  psFree(ids);
2181  //
2182  // Make sure that peaks within each footprint are sorted and unique
2183  //
2184  for (int i = 0; i < footprints->n; ++i) {
2185  pmFootprint *fp = footprints->data[i];
2186  fp->peaks = psArraySort(fp->peaks, pmPeakSortBySN);
2187 
2188  for (int j = 1; j < fp->peaks->n; ++j) { // check for duplicates
2189  if (fp->peaks->data[j] == fp->peaks->data[j-1]) {
2190  (void)psArrayRemoveIndex(fp->peaks, j);
2191  j--; // we moved everything down one
2192  }
2193  }
2194  }
2195 
2196  return PS_ERR_NONE;
2197 }
2198 
2199 /************************************************************************************************************/
2200  /*
2201  * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
2202  * isolated. More precisely, for each peak find the highest coll that you'd have to traverse
2203  * to reach a still higher peak --- and if that coll's more than nsigma DN below your
2204  * starting point, discard the peak.
2205  */
2206 psErrorCode pmFootprintCullPeaks(psImage const *img, // the image wherein lives the footprint
2207  psImage const *weight, // corresponding variance image
2208  pmFootprint *fp, // Footprint containing mortal peaks
2209  float const nsigma_delta, // how many sigma above local background a peak
2210  // needs to be to survive
2211  float const min_threshold) { // minimum permitted coll height
2212  assert (img != NULL);
2213  assert (img->type.type == PS_TYPE_F32);
2214  assert (weight != NULL);
2215  assert (weight->type.type == PS_TYPE_F32);
2216  assert (img->y0 == weight->y0 && img->x0 == weight->x0);
2217  assert (fp != NULL);
2218 
2219  if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
2220  return PS_ERR_NONE;
2221  }
2222 
2223  psRegion subRegion; // desired subregion; 1 larger than bounding box (grr)
2224  subRegion.x0 = fp->bbox.x0;
2225  subRegion.x1 = fp->bbox.x1 + 1;
2226  subRegion.y0 = fp->bbox.y0;
2227  subRegion.y1 = fp->bbox.y1 + 1;
2228  psImage const *subImg = psImageSubset((psImage *)img, subRegion);
2229  psImage const *subWt = psImageSubset((psImage *)weight, subRegion);
2230  assert (subImg != NULL && subWt != NULL);
2231  //
2232  // We need a psArray of peaks brighter than the current peak. We'll fake this
2233  // by reusing the fp->peaks but lying about n.
2234  //
2235  // We do this for efficiency (otherwise I'd need two peaks lists), and we are
2236  // rather too chummy with psArray in consequence. But it works.
2237  //
2238  psArray *brightPeaks = psArrayAlloc(0);
2239  psFree(brightPeaks->data);
2240  brightPeaks->data = psMemIncrRefCounter(fp->peaks->data);// use the data from fp->peaks
2241  //
2242  // The brightest peak is always safe; go through other peaks trying to cull them
2243  //
2244  for (int i = 1; i < fp->peaks->n; ++i) { // n.b. fp->peaks->n can change within the loop
2245  pmPeak const *peak = fp->peaks->data[i];
2246  int x = peak->x - subImg->x0;
2247  int y = peak->y - subImg->y0;
2248  //
2249  // Find the level nsigma below the peak that must separate the peak
2250  // from any of its friends
2251  //
2252  assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
2253  float const stdev = std::sqrt(subWt->data.F32[y][x]);
2254  float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;
2255  if (lsst::utils::isnan(threshold) || threshold < min_threshold) {
2256 #if 1 // min_threshold is assumed to be below the detection threshold,
2257  // so all the peaks are pmFootprint, and this isn't the brightest
2258  (void)psArrayRemoveIndex(fp->peaks, i);
2259  i--; // we moved everything down one
2260  continue;
2261 #else
2262 #error n.b. We will be running LOTS of checks at this threshold, so only find the footprint once
2263  threshold = min_threshold;
2264 #endif
2265  }
2266  if (threshold > subImg->data.F32[y][x]) {
2267  threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
2268  }
2269 
2270  int const peak_id = 1; // the ID for the peak of interest
2271  brightPeaks->n = i; // only stop at a peak brighter than we are
2272  pmFootprint *peakFootprint = pmFindFootprintAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
2273  brightPeaks->n = 0; // don't double free
2274  psImage *idImg = pmSetFootprintID(peakFootprint, peak_id);
2275  psFree(peakFootprint);
2276 
2277  int j;
2278  for (j = 0; j < i; ++j) {
2279  pmPeak const *peak2 = fp->peaks->data[j];
2280  int x2 = peak2->x - subImg->x0;
2281  int y2 = peak2->y - subImg->y0;
2282  int const peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak
2283 
2284  if (peak2_id == peak_id) { // There's a brighter peak within the footprint above
2285  ; // threshold; so cull our initial peak
2286  (void)psArrayRemoveIndex(fp->peaks, i);
2287  i--; // we moved everything down one
2288  break;
2289  }
2290  }
2291  if (j == i) {
2292  ++j;
2293  }
2294 
2295  psFree(idImg);
2296  }
2297 
2298  brightPeaks->n = 0;
2299  psFree(brightPeaks);
2300  psFree((psImage *)subImg);
2301  psFree((psImage *)subWt);
2302 
2303  return PS_ERR_NONE;
2304 }
2305 
2306 /*
2307  * Cull an entire psArray of pmFootprints
2308  */
2309 psErrorCode
2310 pmFootprintArrayCullPeaks(psImage const *img, // the image wherein lives the footprint
2311  psImage const *weight, // corresponding variance image
2312  psArray *footprints, // array of pmFootprints
2313  float const nsigma_delta, // how many sigma above local background a peak
2314  // needs to be to survive
2315  float const min_threshold) { // minimum permitted coll height
2316  for (int i = 0; i < footprints->n; ++i) {
2317  pmFootprint *fp = footprints->data[i];
2318  if (pmFootprintCullPeaks(img, weight, fp, nsigma_delta, min_threshold) != PS_ERR_NONE) {
2319  return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id);
2320  }
2321  }
2322 
2323  return PS_ERR_NONE;
2324 }
2325 
2326 /************************************************************************************************************/
2327 /*
2328  * Extract the peaks in a psArray of pmFootprints, returning a psArray of pmPeaks
2329  */
2330 psArray *pmFootprintArrayToPeaks(psArray const *footprints) {
2331  assert(footprints != NULL);
2332  assert(footprints->n == 0 || pmIsFootprint(footprints->data[0]));
2333 
2334  int npeak = 0;
2335  for (int i = 0; i < footprints->n; ++i) {
2336  pmFootprint const *fp = footprints->data[i];
2337  npeak += fp->peaks->n;
2338  }
2339 
2340  psArray *peaks = psArrayAllocEmpty(npeak);
2341 
2342  for (int i = 0; i < footprints->n; ++i) {
2343  pmFootprint const *fp = footprints->data[i];
2344  for (int j = 0; j < fp->peaks->n; ++j) {
2345  psArrayAdd(peaks, 1, fp->peaks->data[j]);
2346  }
2347  }
2348 
2349  return peaks;
2350 }
2351 #endif
2352 
2353 /************************************************************************************************************/
2354 //
2355 // Explicit instantiations
2356 // \cond
2357 //
2358 //
2359 template
2360 void Footprint::intersectMask(
2361  image::Mask<image::MaskPixel> const& mask,
2362  image::MaskPixel bitMask);
2363 
2364 template
2365 PTR(Footprint) footprintAndMask(
2366  PTR(Footprint) const& foot,
2367  image::Mask<image::MaskPixel>::Ptr const& mask,
2368  image::MaskPixel bitMask);
2369 
2370 template
2372  image::Mask<image::MaskPixel> *mask,
2373  CONST_PTR(std::vector<PTR(Footprint)>) const& footprints,
2374  image::MaskPixel const bitmask);
2375 template
2377  image::Mask<image::MaskPixel> *mask,
2378  std::vector<PTR(Footprint)> const& footprints,
2379  image::MaskPixel const bitmask);
2380 template image::MaskPixel setMaskFromFootprint(
2381  image::Mask<image::MaskPixel> *mask,
2382  Footprint const& foot, image::MaskPixel const bitmask);
2383 template image::MaskPixel clearMaskFromFootprint(
2384  image::Mask<image::MaskPixel> *mask,
2385  Footprint const& foot, image::MaskPixel const bitmask);
2386 
2387 #define INSTANTIATE_NUMERIC(TYPE) \
2388 template \
2389 TYPE setImageFromFootprint(image::Image<TYPE> *image, \
2390  Footprint const& footprint, \
2391  TYPE const value); \
2392 template \
2393 TYPE setImageFromFootprintList(image::Image<TYPE> *image, \
2394  std::vector<PTR(Footprint)> const& footprints, \
2395  TYPE const value); \
2396 template \
2397 TYPE setImageFromFootprintList(image::Image<TYPE> *image, \
2398  CONST_PTR(std::vector<PTR(Footprint)>) footprints, \
2399  TYPE const value); \
2400 template \
2401 void copyWithinFootprint(Footprint const&, \
2402  PTR(lsst::afw::image::Image<TYPE>) const, \
2403  PTR(lsst::afw::image::Image<TYPE>)); \
2404 template \
2405 void copyWithinFootprint(Footprint const&, \
2406  PTR(lsst::afw::image::MaskedImage<TYPE>) const, \
2407  PTR(lsst::afw::image::MaskedImage<TYPE>)); \
2408 template \
2409  void Footprint::clipToNonzero(lsst::afw::image::Image<TYPE> const&); \
2410 
2411 
2412 INSTANTIATE_NUMERIC(float);
2413 INSTANTIATE_NUMERIC(double);
2414 // There's no reason these shouldn't have setImageFromFootprint(), etc, instantiated
2415 INSTANTIATE_NUMERIC(boost::uint16_t);
2416 INSTANTIATE_NUMERIC(int);
2417 INSTANTIATE_NUMERIC(boost::uint64_t);
2418 
2419 
2420 #define INSTANTIATE_MASK(PIXEL) \
2421 template \
2422 void Footprint::insertIntoImage( \
2423  lsst::afw::image::Image<PIXEL>& idImage, \
2424  boost::uint64_t const id, \
2425  geom::Box2I const& region=geom::Box2I() \
2426  ) const; \
2427 template \
2428 void Footprint::insertIntoImage( \
2429  lsst::afw::image::Image<PIXEL>& idImage, \
2430  boost::uint64_t const id, \
2431  bool const overwriteId, long const idMask, \
2432  std::set<boost::uint64_t> *oldIds, \
2433  geom::Box2I const& region=geom::Box2I() \
2434  ) const; \
2435 
2436 INSTANTIATE_MASK(boost::uint16_t);
2437 INSTANTIATE_MASK(int);
2438 INSTANTIATE_MASK(boost::uint64_t);
2439 
2440 
2441 }}}
2442 // \endcond
2443 
2444 // LocalWords: SpanList
int y
int iter
void nearestFootprint(std::vector< boost::shared_ptr< Footprint >> const &foots, lsst::afw::image::Image< boost::uint16_t >::Ptr argmin, lsst::afw::image::Image< boost::uint16_t >::Ptr dist)
boost::uint16_t MaskPixel
Declare the Kernel class and subclasses.
A coordinate class intended to represent absolute positions.
#define PTR(...)
Definition: base.h:41
Represent a set of pixels of an arbitrary shape and size.
table::Key< std::string > name
Definition: ApCorrMap.cc:71
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
tbl::Key< double > weight
void swap(ImageBase< PixelT > &a, ImageBase< PixelT > &b)
Definition: Image.cc:291
void reserve(size_type n)
Increase the capacity of the catalog to the given size.
Definition: Catalog.h:421
CatalogT< BaseRecord > BaseCatalog
Definition: fwd.h:61
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:47
Box2I BoxI
Definition: Box.h:479
geom::Point2D skyToPixel(geom::Angle sky1, geom::Angle sky2) const
Convert from sky coordinates (e.g. RA/dec) to pixel positions.
Definition: Wcs.cc:782
boost::shared_ptr< Footprint > footprintAndMask(boost::shared_ptr< Footprint > const &foot, typename image::Mask< MaskT >::Ptr const &mask, MaskT const bitmask)
Return a Footprint that&#39;s the intersection of a Footprint with a Mask.
#define CONST_PTR(...)
Definition: base.h:47
lsst::afw::detection::Footprint Footprint
Definition: Source.h:61
definition of the Trace messaging facilities
float _value
Definition: saturated.cc:31
SelectEigenView< T >::Type copy(Eigen::EigenBase< T > const &other)
Copy an arbitrary Eigen expression into a new EigenView.
Definition: eigen.h:390
void include(Point2D const &point)
Expand this to ensure that this-&gt;contains(point).
int getMinY() const
Definition: Box.h:125
Implementation of the WCS standard for a any projection.
Definition: Wcs.h:107
void shift(Vector< T, N > const &offset, Array< std::complex< T >, N, C > const &array, int const real_last_dim)
Perform a Fourier-space translation transform.
Definition: FourierOps.h:146
afw::table::CatalogT< PeakRecord > PeakCatalog
Definition: Peak.h:225
Point< int, 2 > Point2I
Definition: PSF.h:39
boost::shared_ptr< coord::Coord > pixelToSky(double pix1, double pix2) const
Convert from pixel position to sky coordinates (e.g. RA/dec)
Definition: Wcs.cc:858
int getMaxX() const
Definition: Box.h:128
int isnan(T t)
Definition: ieee.h:110
int const x0
Definition: saturated.cc:45
An integer coordinate rectangle.
Definition: Box.h:53
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
Represent a collections of footprints associated with image data.
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition: BaseRecord.h:123
boost::shared_ptr< Footprint > shrinkFootprint(Footprint const &foot, int nGrow, bool isotropic)
afw::geom::ellipses::Quadrupole Shape
Definition: constants.h:58
void shift(Extent2I const &offset)
Shift the position of the box by the given offset.
int getWidth() const
Return the number of columns in the image.
Definition: Image.h:237
int getHeight() const
Definition: Box.h:155
boost::shared_ptr< Footprint > growFootprint(Footprint const &foot, int nGrow, bool isotropic=true)
CatalogIterator< typename Internal::const_iterator > const_iterator
Definition: Catalog.h:108
Extent2I const getDimensions() const
Definition: Box.h:153
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:93
geom::Box2I const & _bbox
Definition: fits_io_mpl.h:80
x_iterator x_at(int x, int y) const
Return an x_iterator to the point (x, y) in the image.
Definition: Image.h:329
if(width!=gim.getWidth()||height!=gim.getHeight()||x0!=gim.getX0()||y0!=gim.getY0())
Definition: saturated.cc:47
Footprint(int nspan=0, geom::Box2I const &region=geom::Box2I())
bool contains(Point2I const &point) const
Return true if the box contains the point.
int getMaxY() const
Definition: Box.h:129
int getMinX() const
Definition: Box.h:124
void include(Point2I const &point)
Expand this to ensure that this-&gt;contains(point).
int getY0() const
Definition: Image.h:255
MaskT clearMaskFromFootprint(lsst::afw::image::Mask< MaskT > *mask, Footprint const &footprint, MaskT const bitmask)
(AND ~bitmask) all the Mask&#39;s pixels that are in the Footprint; that is, set to zero in the Mask-inte...
LSST bitmasks.
int x
ImageT::Pixel setImageFromFootprint(ImageT *image, Footprint const &footprint, typename ImageT::Pixel const value)
Set all image pixels in a Footprint to a given value.
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
std::vector< lsst::afw::geom::Box2I > footprintToBBoxList(Footprint const &foot)
bool isEmpty() const
Return true if the box contains no points.
Definition: Box.h:166
tuple m
Definition: lsstimport.py:48
Extent< int, N > floor(Extent< double, N > const &input)
geom::Box2I getBBox(ImageOrigin origin=PARENT) const
Definition: Image.h:377
MaskT setMaskFromFootprint(lsst::afw::image::Mask< MaskT > *mask, Footprint const &footprint, MaskT const bitmask)
OR bitmask into all the Mask&#39;s pixels that are in the Footprint.
int id
Definition: CR.cc:151
afw::table::Key< double > b
Point< double, 2 > Point2D
Definition: Point.h:286
xy_locator xy_at(int x, int y) const
Definition: Image.h:351
int getHeight() const
Return the number of rows in the image.
Definition: Image.h:239
A floating-point coordinate rectangle geometry.
Definition: Box.h:271
boost::shared_ptr< Footprint > mergeFootprints(Footprint const &foot1, Footprint const &foot2)
int const y0
Definition: saturated.cc:45
Record class that represents a peak in a Footprint.
Definition: Peak.h:40
ImageT val
Definition: CR.cc:154
Extent< int, 2 > Extent2I
Definition: Extent.h:355
int getWidth() const
Definition: Box.h:154
Utility functions for kernels.
void copyWithinFootprint(Footprint const &foot, boost::shared_ptr< ImageOrMaskedImageT > const input, boost::shared_ptr< ImageOrMaskedImageT > output)
ImageT::Pixel setImageFromFootprintList(ImageT *image, boost::shared_ptr< std::vector< boost::shared_ptr< Footprint >> const > footprints, typename ImageT::Pixel const value)
Set all image pixels in a set of Footprints to a given value.
MaskT setMaskFromFootprintList(lsst::afw::image::Mask< MaskT > *mask, std::vector< boost::shared_ptr< Footprint >> const &footprints, MaskT const bitmask)
OR bitmask into all the Mask&#39;s pixels which are in the set of Footprints.
Include files required for standard LSST Exception handling.
int getX0() const
Definition: Image.h:247