LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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"
38
39namespace lsst {
40namespace meas {
41namespace algorithms {
42
43/************************************************************************************************************/
44/*
45 * PsfCandidate's members
46 */
47template <typename PixelT>
48int PsfCandidate<PixelT>::_border = 0;
49template <typename PixelT>
50int PsfCandidate<PixelT>::_defaultWidth = 21;
51template <typename PixelT>
52float PsfCandidate<PixelT>::_pixelThreshold = 0.0;
53template <typename PixelT>
54bool PsfCandidate<PixelT>::_doMaskBlends = true;
55
56/************************************************************************************************************/
57namespace {
58template <typename T> // functor used by makeImageFromMask to return inputMask
59struct noop : public afw::image::pixelOp1<T> {
60 T operator()(T x) const { return x; }
61};
62
63template <typename T> // functor used by makeImageFromMask to return (inputMask & mask)
64struct andMask : public afw::image::pixelOp1<T> {
65 andMask(T mask) : _mask(mask) {}
66 T operator()(T x) const { return (x & _mask); }
67
68private:
69 T _mask;
70};
71
72template <typename T>
73andMask<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 */
80template <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
103double 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
111template <typename MaskT>
112class BlendedFunctor {
113public:
114 BlendedFunctor(afw::detection::PeakRecord const& central,
115 afw::detection::PeakCatalog const& peaks,
116 afw::image::MaskPixel turnOff,
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
136private:
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
163template <typename PixelT>
165PsfCandidate<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
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 std::shared_ptr<afw::detection::Footprint const> foot = getSource()->getFootprint();
194 typedef afw::detection::PeakCatalog PeakCatalog;
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
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) {
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 */
228
229 std::shared_ptr<afw::image::Image<int>> mim = makeImageFromMask<int>(*image->getMask(), makeAndMask(detected));
231 fs = std::make_shared<afw::detection::FootprintSet>(*mim, afw::detection::Threshold(1));
232 std::shared_ptr<FootprintList const> 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) {
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) {
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) {
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
275template <typename PixelT>
277PsfCandidate<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
289template <typename PixelT>
292 int const width = getWidth() == 0 ? _defaultWidth : getWidth();
293 int const height = getHeight() == 0 ? _defaultWidth : getHeight();
294 return getMaskedImage(width, height);
295}
296
297template <typename PixelT>
299 return _border;
300}
301
302template <typename PixelT>
304 _border = border;
305}
306
307template <typename PixelT>
309 _pixelThreshold = threshold;
310}
311
312template <typename PixelT>
314 return _pixelThreshold;
315}
316
317template <typename PixelT>
318void PsfCandidate<PixelT>::setMaskBlends(bool doMaskBlends) {
319 _doMaskBlends = doMaskBlends;
320}
321
322template <typename PixelT>
324 return _doMaskBlends;
325}
326
332template <typename PixelT>
334PsfCandidate<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 std::shared_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 std::shared_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//
364typedef float Pixel;
365// template class PsfCandidate<afw::image::MaskedImage<Pixel> >;
366template class PsfCandidate<Pixel>;
368
369} // namespace algorithms
370} // namespace meas
371} // namespace lsst
AmpInfoBoxKey bbox
Definition Amplifier.cc:117
int end
#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:95
int y
Definition SpanSet.cc:48
std::vector< std::shared_ptr< Footprint > > FootprintList
The FootprintSet's set of Footprints.
@ PIXEL_STDEV
Use number of sigma given per-pixel s.d.
Definition Threshold.h:51
typename _const_view_t::x_iterator const_x_iterator
A const iterator for traversing the pixels in a row.
Definition ImageBase.h:141
typename _view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition ImageBase.h:133
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.
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.
std::shared_ptr< afw::image::MaskedImage< PixelT > const > getMaskedImage() const
Return the image at the position of the Source, without any sub-pixel shifts to put the centre of the...
static void setMaskBlends(bool doMaskBlends)
Set whether blends are masked.
std::shared_ptr< afw::image::MaskedImage< PixelT > > getOffsetImage(std::string const algorithm, unsigned int buffer) const
Return an offset version of the image of the source.
T infinity(T... args)
afw::table::CatalogT< PeakRecord > PeakCatalog
Definition Peak.h:244
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.
T pow(T... args)
ImageT val
Definition CR.cc:146