LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
PsfCandidate.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008-2015 AURA/LSST.
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 
33 #include "lsst/geom.h"
35 #include "lsst/afw/image/Image.h"
38 
39 namespace lsst {
40 namespace meas {
41 namespace algorithms {
42 
43 /************************************************************************************************************/
44 /*
45  * PsfCandidate's members
46  */
47 template <typename PixelT>
48 int PsfCandidate<PixelT>::_border = 0;
49 template <typename PixelT>
50 int PsfCandidate<PixelT>::_defaultWidth = 21;
51 template <typename PixelT>
52 float PsfCandidate<PixelT>::_pixelThreshold = 0.0;
53 template <typename PixelT>
54 bool PsfCandidate<PixelT>::_doMaskBlends = true;
55 
56 /************************************************************************************************************/
57 namespace {
58 template <typename T> // functor used by makeImageFromMask to return inputMask
59 struct noop : public afw::image::pixelOp1<T> {
60  T operator()(T x) const { return x; }
61 };
62 
63 template <typename T> // functor used by makeImageFromMask to return (inputMask & mask)
64 struct andMask : public afw::image::pixelOp1<T> {
65  andMask(T mask) : _mask(mask) {}
66  T operator()(T x) const { return (x & _mask); }
67 
68 private:
69  T _mask;
70 };
71 
72 template <typename T>
73 andMask<T> makeAndMask(T val) {
74  return andMask<T>(val);
75 }
76 
77 /*
78  * Return an Image initialized from a Mask (possibly modified by func)
79  */
80 template <typename LhsT, typename RhsT>
82  afw::image::Mask<RhsT> const& rhs,
83  afw::image::pixelOp1<RhsT> const& func = noop<RhsT>()
84 ) {
86  std::make_shared<afw::image::Image<LhsT>>(rhs.getDimensions());
87  lhs->setXY0(rhs.getXY0());
88 
89  for (int y = 0; y != lhs->getHeight(); ++y) {
90  typename afw::image::Image<RhsT>::const_x_iterator rhsPtr = rhs.row_begin(y);
91 
92  for (typename afw::image::Image<LhsT>::x_iterator lhsPtr = lhs->row_begin(y),
93  lhsEnd = lhs->row_end(y);
94  lhsPtr != lhsEnd; ++rhsPtr, ++lhsPtr) {
95  *lhsPtr = func(*rhsPtr);
96  }
97  }
98 
99  return lhs;
100 }
101 
103 double distanceSquared(double x, double y, afw::detection::PeakRecord const& peak) {
104  return std::pow(peak.getIx() - x, 2) + std::pow(peak.getIy() - y, 2);
105 }
106 
111 template <typename MaskT>
112 class BlendedFunctor {
113 public:
114  BlendedFunctor(afw::detection::PeakRecord const& central,
115  afw::detection::PeakCatalog const& peaks,
116  afw::image::MaskPixel turnOff,
117  afw::image::MaskPixel turnOn
118  )
119  : _central(central), _peaks(peaks), _turnOff(~turnOff), _turnOn(turnOn) {}
120 
122  void operator()(geom::Point2I const& point, MaskT& val) {
123  int x = point.getX();
124  int y = point.getY();
125  double const central = distanceSquared(x, y, _central);
126  for (afw::detection::PeakCatalog::const_iterator iter = _peaks.begin(), end = _peaks.end();
127  iter != end; ++iter) {
128  double const dist2 = distanceSquared(x, y, *iter);
129  if (dist2 < central) {
130  val &= _turnOff;
131  val |= _turnOn;
132  }
133  }
134  }
135 
136 private:
137  afw::detection::PeakRecord const& _central;
138  afw::detection::PeakCatalog const& _peaks;
139  afw::image::MaskPixel const _turnOff;
140  afw::image::MaskPixel const _turnOn;
141 };
142 
143 } // anonymous namespace
144 
163 template <typename PixelT>
164 PTR(afw::image::MaskedImage<PixelT>)
165 PsfCandidate<PixelT>::extractImage(unsigned int width, // Width of image
166  unsigned int height // Height of image
167  ) const {
168  geom::Point2I const cen(afw::image::positionToIndex(getXCenter()),
169  afw::image::positionToIndex(getYCenter()));
170  geom::Point2I const llc(cen[0] - width / 2 - _parentExposure->getX0(),
171  cen[1] - height / 2 - _parentExposure->getY0());
172 
173  geom::BoxI bbox(llc, geom::ExtentI(width, height));
174 
175  PTR(MaskedImageT) image;
176  try {
177  MaskedImageT mimg = _parentExposure->getMaskedImage();
178  image.reset(new MaskedImageT(mimg, bbox, afw::image::LOCAL, true)); // a deep copy
179  } catch (pex::exceptions::LengthError& e) {
180  LSST_EXCEPT_ADD(e, "Extracting image of PSF candidate");
181  throw e;
182  }
183 
184  //
185  // Set INTRP and unset DETECTED for any pixels we don't want to deal with.
186  //
187  afw::image::MaskPixel const intrp =
188  MaskedImageT::Mask::getPlaneBitMask("INTRP"); // mask bit for bad pixels
189  afw::image::MaskPixel const detected = MaskedImageT::Mask::getPlaneBitMask("DETECTED"); // object pixels
190 
191  // Mask out blended objects
192  if (getMaskBlends()) {
193  CONST_PTR(afw::detection::Footprint) foot = getSource()->getFootprint();
195  PeakCatalog const& peaks = foot->getPeaks();
196  if (peaks.size() > 1) {
197  // Mask all pixels in the footprint except for those closest to the central peak
199  PTR(afw::detection::PeakRecord) central;
200  for (PeakCatalog::const_iterator iter = peaks.begin(), end = peaks.end(); iter != end; ++iter) {
201  double const dist2 = distanceSquared(getXCenter(), getYCenter(), *iter);
202  if (dist2 < best) {
203  best = dist2;
204  central = iter;
205  }
206  }
207  assert(central); // We must have found something
208 
209  PeakCatalog others(peaks.getTable());
210  others.reserve(peaks.size() - 1);
211  for (PeakCatalog::const_iterator iter = peaks.begin(), end = peaks.end(); iter != end; ++iter) {
212  PTR(afw::detection::PeakRecord) ptr(iter);
213  if (central != ptr) {
214  others.push_back(ptr);
215  }
216  }
217 
218  BlendedFunctor<typename MaskedImageT::Mask::Pixel> functor(*central, others, detected, intrp);
219  foot->getSpans()->clippedTo(image->getBBox())->applyFunctor(functor, *image->getMask());
220  }
221  }
222 
223  /*
224  * Mask any DETECTED pixels other than the one in the center of the object;
225  * we grow the Footprint a bit first
226  */
227  typedef afw::detection::FootprintSet::FootprintList FootprintList;
228 
229  PTR(afw::image::Image<int>) mim = makeImageFromMask<int>(*image->getMask(), makeAndMask(detected));
230  PTR(afw::detection::FootprintSet)
231  fs = std::make_shared<afw::detection::FootprintSet>(*mim, afw::detection::Threshold(1));
232  CONST_PTR(FootprintList) feet = fs->getFootprints();
233 
234  if (feet->size() > 1) {
235  int const ngrow = 3; // number of pixels to grow bad Footprints
236  //
237  // Go through Footprints looking for ones that don't contain cen
238  //
239  for (FootprintList::const_iterator fiter = feet->begin(); fiter != feet->end(); ++fiter) {
240  PTR(afw::detection::Footprint) foot = *fiter;
241  if (foot->contains(cen)) {
242  continue;
243  }
244 
245  // Dilate and clip to the image bounding box, incase the span grows outside the image
246  auto bigSpan = foot->getSpans()->dilated(ngrow)->clippedTo(image->getBBox());
247  bigSpan->clearMask(*image->getMask(), detected);
248  bigSpan->setMask(*image->getMask(), intrp);
249  }
250  }
251 
252  // Mask high pixels unconnected to the center
253  if (_pixelThreshold > 0.0) {
254  CONST_PTR(afw::detection::FootprintSet)
255  fpSet = std::make_shared<afw::detection::FootprintSet>(
256  *image, afw::detection::Threshold(_pixelThreshold, afw::detection::Threshold::PIXEL_STDEV));
257  for (FootprintList::const_iterator fpIter = fpSet->getFootprints()->begin();
258  fpIter != fpSet->getFootprints()->end(); ++fpIter) {
259  CONST_PTR(afw::detection::Footprint) fp = *fpIter;
260  if (!fp->contains(cen)) {
261  fp->getSpans()->clearMask(*image->getMask(), detected);
262  fp->getSpans()->setMask(*image->getMask(), intrp);
263  }
264  }
265  }
266 
267  return image;
268 }
269 
275 template <typename PixelT>
276 CONST_PTR(afw::image::MaskedImage<PixelT>)
277 PsfCandidate<PixelT>::getMaskedImage(int width, int height) const {
278  if (!_image || (width != _image->getWidth() || height != _image->getHeight())) {
279  _image = extractImage(width, height);
280  }
281  return _image;
282 }
283 
289 template <typename PixelT>
291 PsfCandidate<PixelT>::getMaskedImage() const {
292  int const width = getWidth() == 0 ? _defaultWidth : getWidth();
293  int const height = getHeight() == 0 ? _defaultWidth : getHeight();
294  return getMaskedImage(width, height);
295 }
296 
297 template <typename PixelT>
299  return _border;
300 }
301 
302 template <typename PixelT>
304  _border = border;
305 }
306 
307 template <typename PixelT>
309  _pixelThreshold = threshold;
310 }
311 
312 template <typename PixelT>
314  return _pixelThreshold;
315 }
316 
317 template <typename PixelT>
318 void PsfCandidate<PixelT>::setMaskBlends(bool doMaskBlends) {
319  _doMaskBlends = doMaskBlends;
320 }
321 
322 template <typename PixelT>
324  return _doMaskBlends;
325 }
326 
332 template <typename PixelT>
334 PsfCandidate<PixelT>::getOffsetImage(std::string const algorithm, // Warping algorithm to use
335  unsigned int buffer // Buffer for warping
336  ) const {
337  unsigned int const width = getWidth() == 0 ? _defaultWidth : getWidth();
338  unsigned int const height = getHeight() == 0 ? _defaultWidth : getHeight();
339  if (_offsetImage && static_cast<unsigned int>(_offsetImage->getWidth()) == width + 2 * buffer &&
340  static_cast<unsigned int>(_offsetImage->getHeight()) == height + 2 * buffer) {
341  return _offsetImage;
342  }
343 
344  PTR(MaskedImageT) image = extractImage(width + 2 * buffer, height + 2 * buffer);
345 
346  double const xcen = getXCenter(), ycen = getYCenter();
347  double const dx = afw::image::positionToIndex(xcen, true).second;
348  double const dy = afw::image::positionToIndex(ycen, true).second;
349 
350  PTR(MaskedImageT) offset = afw::math::offsetImage(*image, -dx, -dy, algorithm);
351  geom::Point2I llc(buffer, buffer);
352  geom::Extent2I dims(width, height);
353  geom::Box2I box(llc, dims);
354  _offsetImage.reset(new MaskedImageT(*offset, box, afw::image::LOCAL, true)); // Deep copy
355 
356  return _offsetImage;
357 }
358 
359 /************************************************************************************************************/
360 //
361 // Explicit instantiations
362 //
364 typedef float Pixel;
365 // template class PsfCandidate<afw::image::MaskedImage<Pixel> >;
366 template class PsfCandidate<Pixel>;
368 
369 } // namespace algorithms
370 } // namespace meas
371 } // namespace lsst
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
int end
double x
#define LSST_EXCEPT_ADD(e, m)
Add the current location and a message to an existing exception before rethrowing it.
Definition: Exception.h:54
afw::table::Key< afw::table::Array< ImagePixelT > > image
afw::table::Key< afw::table::Array< MaskPixelT > > mask
Class used by SpatialCell for spatial PSF fittig.
uint64_t * ptr
Definition: RangeSet.cc:88
int y
Definition: SpanSet.cc:49
#define PTR(...)
Definition: base.h:41
#define CONST_PTR(...)
A shared pointer to a const object.
Definition: base.h:47
std::vector< std::shared_ptr< Footprint > > FootprintList
The FootprintSet's set of Footprints.
Definition: FootprintSet.h:56
@ PIXEL_STDEV
Use number of sigma given per-pixel s.d.
Definition: Threshold.h:51
_const_view_t::x_iterator const_x_iterator
A const iterator for traversing the pixels in a row.
Definition: ImageBase.h:141
_view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: ImageBase.h:133
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
CatalogIterator< typename Internal::const_iterator > const_iterator
Definition: Catalog.h:111
An integer coordinate rectangle.
Definition: Box.h:55
Class stored in SpatialCells for spatial Psf fitting.
Definition: PsfCandidate.h:55
static float getPixelThreshold()
Get threshold for rejecting pixels unconnected with the central footprint.
static bool getMaskBlends()
Get whether blends are masked.
static void setBorderWidth(int border)
Set the number of pixels to ignore around the candidate image's edge.
static int getBorderWidth()
Return the number of pixels being ignored around the candidate image's edge.
static void setPixelThreshold(float threshold)
Set threshold for rejecting pixels unconnected with the central footprint.
static void setMaskBlends(bool doMaskBlends)
Set whether blends are masked.
T infinity(T... args)
afw::table::CatalogT< PeakRecord > PeakCatalog
Definition: Peak.h:244
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
std::int32_t MaskPixel
default type for Masks and MaskedImage Masks
int positionToIndex(double pos)
Convert image position to nearest integer index.
Definition: ImageUtils.h:69
std::shared_ptr< ImageT > offsetImage(ImageT const &image, float dx, float dy, std::string const &algorithmName="lanczos5", unsigned int buffer=0)
Return an image offset by (dx, dy) using the specified algorithm.
Definition: offsetImage.cc:41
lsst::afw::detection::Footprint Footprint
Definition: Source.h:59
float Pixel
Typedefs to be used for pixel values.
Definition: common.h:37
A base class for image defects.
STL namespace.
T pow(T... args)
ImageT val
Definition: CR.cc:146