LSST Applications 26.0.0,g0265f82a02+6660c170cc,g07994bdeae+30b05a742e,g0a0026dc87+17526d298f,g0a60f58ba1+17526d298f,g0e4bf8285c+96dd2c2ea9,g0ecae5effc+c266a536c8,g1e7d6db67d+6f7cb1f4bb,g26482f50c6+6346c0633c,g2bbee38e9b+6660c170cc,g2cc88a2952+0a4e78cd49,g3273194fdb+f6908454ef,g337abbeb29+6660c170cc,g337c41fc51+9a8f8f0815,g37c6e7c3d5+7bbafe9d37,g44018dc512+6660c170cc,g4a941329ef+4f7594a38e,g4c90b7bd52+5145c320d2,g58be5f913a+bea990ba40,g635b316a6c+8d6b3a3e56,g67924a670a+bfead8c487,g6ae5381d9b+81bc2a20b4,g93c4d6e787+26b17396bd,g98cecbdb62+ed2cb6d659,g98ffbb4407+81bc2a20b4,g9ddcbc5298+7f7571301f,ga1e77700b3+99e9273977,gae46bcf261+6660c170cc,gb2715bf1a1+17526d298f,gc86a011abf+17526d298f,gcf0d15dbbd+96dd2c2ea9,gdaeeff99f8+0d8dbea60f,gdb4ec4c597+6660c170cc,ge23793e450+96dd2c2ea9,gf041782ebf+171108ac67
LSST Data Management Base Package
Loading...
Searching...
No Matches
Blendedness.cc
Go to the documentation of this file.
1// -*- LSST-C++ -*-
2/*
3 * This file is part of meas_base.
4 *
5 * Developed for the LSST Data Management System.
6 * This product includes software developed by the LSST Project
7 * (https://www.lsst.org).
8 * See the COPYRIGHT file at the top-level directory of this distribution
9 * for details of code ownership.
10 *
11 * This program is free software: you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License as published by
13 * the Free Software Foundation, either version 3 of the License, or
14 * (at your option) any later version.
15 *
16 * This program is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with this program. If not, see <https://www.gnu.org/licenses/>.
23 */
24
25#include <cmath>
26
27#include "boost/math/constants/constants.hpp"
28
35
36namespace lsst {
37namespace meas {
38namespace base {
39namespace {
40FlagDefinitionList flagDefinitions;
41} // namespace
42
43FlagDefinition const BlendednessAlgorithm::FAILURE = flagDefinitions.addFailureFlag();
44FlagDefinition const BlendednessAlgorithm::NO_CENTROID =
45 flagDefinitions.add("flag_noCentroid", "Object has no centroid");
46FlagDefinition const BlendednessAlgorithm::NO_SHAPE =
47 flagDefinitions.add("flag_noShape", "Object has no shape");
48
50
51namespace {
52
53double computeOldBlendedness(std::shared_ptr<afw::detection::Footprint const> childFootprint,
54 afw::image::Image<float> const& parentImage) {
55 if (!childFootprint) {
56 throw LSST_EXCEPT(pex::exceptions::LogicError, "blendedness_old requires a Footprint.");
57 }
58
60 childHeavy = std::dynamic_pointer_cast<afw::detection::HeavyFootprint<float> const>(childFootprint);
61
62 if (!childHeavy) {
63 return 0.0; // if it's not a HeavyFootprint, it's not blended.
64 }
65
66 if (!parentImage.getBBox(afw::image::PARENT).contains(childHeavy->getBBox())) {
67 throw LSST_EXCEPT(pex::exceptions::LogicError, "Child footprint extends beyond image.");
68 }
69
70 // Iterate over all the spans in the child HeavyFootprint,
71 // along with iterators for the child pixels (from the HeavyFootprint)
72 // and parent pixels (from the Exposure).
74 typedef ndarray::Array<float const, 1, 1>::Iterator ChildPixIter;
75 auto spanIter = childHeavy->getSpans()->begin();
76 auto const spanEnd = childHeavy->getSpans()->end();
77 ChildPixIter childPixIter = childHeavy->getImageArray().begin();
78 double cp = 0.0; // child.dot(parent)
79 double cc = 0.0; // child.dot(child)
80 while (spanIter != spanEnd) {
81 afw::geom::Span const& span = *spanIter;
82 ParentPixIter parentPixIter =
83 parentImage.x_at(span.getBeginX() - parentImage.getX0(), span.getY() - parentImage.getY0());
84 int const width = span.getWidth();
85 // Iterate over the pixels within the span, updating the dot products.
86 for (int x = 0; x < width; ++parentPixIter, ++childPixIter, ++x) {
87 cp += (*childPixIter) * ((*parentPixIter) - (*childPixIter));
88 cc += (*childPixIter) * (*childPixIter);
89 }
90 ++spanIter;
91 }
92 if (cc > 0.0) {
93 return cp / cc;
94 }
95 return 0.0;
96}
97
98class FluxAccumulator {
99public:
100 FluxAccumulator() : _w(0.0), _ww(0.0), _wd(0.0) {}
101
102 void operator()(double, double, float weight, float data) {
103 _w += weight;
104 _ww += weight * weight;
105 _wd += weight * data;
106 }
107
108 double getFlux() const { return _w * _wd / _ww; }
109
110protected:
111 double _w;
112 double _ww;
113 double _wd;
114};
115
116class ShapeAccumulator : public FluxAccumulator {
117public:
118 ShapeAccumulator() : FluxAccumulator(), _wdxx(0.0), _wdyy(0.0), _wdxy(0.0) {}
119
120 void operator()(double x, double y, float weight, float data) {
121 FluxAccumulator::operator()(x, y, weight, data);
122 _wdxx += x * x * weight * data;
123 _wdyy += y * y * weight * data;
124 _wdxy += x * y * weight * data;
125 }
126
127 ShapeResult getShape() const {
128 // Factor of 2 corrects for bias from weight function (correct is exact for an object
129 // with a Gaussian profile.)
130 ShapeResult result;
131 result.xx = 2.0 * _wdxx / _wd;
132 result.yy = 2.0 * _wdyy / _wd;
133 result.xy = 2.0 * _wdxy / _wd;
134 return result;
135 }
136
137private:
138 double _wdxx;
139 double _wdyy;
140 double _wdxy;
141};
142
143template <typename Accumulator>
144void computeMoments(afw::image::MaskedImage<float> const& image, geom::Point2D const& centroid,
145 afw::geom::ellipses::Quadrupole const& shape, double nSigmaWeightMax,
146 Accumulator& accumulatorRaw, Accumulator& accumulatorAbs) {
147 geom::Box2I bbox = image.getBBox(afw::image::PARENT);
148
149 afw::geom::ellipses::Ellipse ellipse(shape, centroid);
150 ellipse.getCore().scale(nSigmaWeightMax);
151
152 // To evaluate an elliptically-symmetric function, we transform points
153 // by the following transform, then evaluate a circularly-symmetric function
154 // at the transformed positions.
155 geom::LinearTransform transform = shape.getGridTransform();
156
157 typedef afw::geom::ellipses::PixelRegion::Iterator SpanIter; // yields Spans
158 typedef afw::geom::Span::Iterator PointIter; // yields Point2I positions
159 typedef afw::image::MaskedImage<float>::const_x_iterator PixelIter; // yields pixel values
160
161 afw::geom::ellipses::PixelRegion region(ellipse);
162 bool isContained = bbox.contains(region.getBBox());
163 SpanIter const spanEnd = region.end();
164 for (SpanIter spanIter = region.begin(); spanIter != spanEnd; ++spanIter) {
165 afw::geom::Span span = *spanIter;
166 if (!isContained) {
167 if (span.getY() < bbox.getMinY() || span.getY() > bbox.getMaxY()) {
168 continue;
169 }
170 span = afw::geom::Span(span.getY(), std::max(span.getMinX(), bbox.getMinX()),
171 std::min(span.getMaxX(), bbox.getMaxX()));
172 if (span.getMinX() > span.getMaxX()) {
173 continue;
174 }
175 }
176 PixelIter pixelIter = image.x_at(span.getBeginX() - image.getX0(), span.getY() - image.getY0());
177 PointIter const pointEnd = span.end();
178 for (PointIter pointIter = span.begin(); pointIter != pointEnd; ++pointIter, ++pixelIter) {
179 geom::Extent2D d = geom::Point2D(*pointIter) - centroid;
181 // use single precision for faster exp, erf
182 float weight = std::exp(static_cast<float>(-0.5 * td.computeSquaredNorm()));
183 float data = pixelIter.image();
184 accumulatorRaw(d.getX(), d.getY(), weight, data);
185 float variance = pixelIter.variance();
188 accumulatorAbs(d.getX(), d.getY(), weight, std::abs(data) - bias);
189 }
190 }
191}
192
193} // namespace
194
197 : _ctrl(ctrl) {
198 if (_ctrl.doOld) {
199 _old = schema.addField<double>(
200 schema.join(name, "old"),
201 "Blendedness from dot products: (child.dot(parent)/child.dot(child) - 1)");
202 }
203 if (_ctrl.doFlux) {
204 _raw = schema.addField<double>(
205 schema.join(name, "raw"),
206 "Measure of how much the flux is affected by neighbors: "
207 "(1 - child_instFlux/parent_instFlux). Operates on the \"raw\" pixel values.");
208 _instFluxChildRaw = schema.addField<double>(
209 schema.join(name, "raw_child_instFlux"),
210 "Instrumental flux of the child, measured with a Gaussian weight matched to the child. "
211 "Operates on the \"raw\" pixel values.", "count");
212 _instFluxParentRaw = schema.addField<double>(
213 schema.join(name, "raw_parent_instFlux"),
214 "Instrumental flux of the parent, measured with a Gaussian weight matched to the child. "
215 "Operates on the \"raw\" pixel values.", "count");
216 _abs = schema.addField<double>(
217 schema.join(name, "abs"),
218 "Measure of how much the flux is affected by neighbors: "
219 "(1 - child_instFlux/parent_instFlux). "
220 "Operates on the absolute value of the pixels to try to obtain a \"de-noised\" value. "
221 "See section 4.9.11 of Bosch et al. 2018, PASJ, 70, S5 for details.");
222 _instFluxChildAbs = schema.addField<double>(
223 schema.join(name, "abs_child_instFlux"),
224 "Instrumental flux of the child, measured with a Gaussian weight matched to the child. "
225 "Operates on the absolute value of the pixels to try to obtain a \"de-noised\" value. "
226 "See section 4.9.11 of Bosch et al. 2018, PASJ, 70, S5 for details.", "count");
227 _instFluxParentAbs = schema.addField<double>(
228 schema.join(name, "abs_parent_instFlux"),
229 "Instrumental flux of the parent, measured with a Gaussian weight matched to the child. "
230 "Operates on the absolute value of the pixels to try to obtain a \"de-noised\" value. "
231 "See section 4.9.11 of Bosch et al. 2018, PASJ, 70, S5 for details.", "count");
232 }
233 if (_ctrl.doShape) {
234 _shapeChildRaw = ShapeResultKey::addFields(
235 schema, schema.join(name, "raw_child"),
236 "Shape of the child, measured with a Gaussian weight matched to the child. "
237 "Operates on the \"raw\" pixel values.", NO_UNCERTAINTY);
238 _shapeParentRaw = ShapeResultKey::addFields(
239 schema, schema.join(name, "raw_parent"),
240 "Shape of the parent, measured with a Gaussian weight matched to the child. "
241 "Operates on the \"raw\" pixel values.", NO_UNCERTAINTY);
242 _shapeChildAbs = ShapeResultKey::addFields(
243 schema, schema.join(name, "abs_child"),
244 "Shape of the child, measured with a Gaussian weight matched to the child. "
245 "Operates on the absolute value of the pixels to try to obtain a \"de-noised\" value. "
246 "See section 4.9.11 of Bosch et al. 2018, PASJ, 70, S5 for details.",
248 _shapeParentAbs = ShapeResultKey::addFields(
249 schema, schema.join(name, "abs_parent"),
250 "Shape of the parent, measured with a Gaussian weight matched to the child. "
251 "Operates on the absolute value of the pixels to try to obtain a \"de-noised\" value. "
252 "See section 4.9.11 of Bosch et al. 2018, PASJ, 70, S5 for details.",
254 }
255 if (_ctrl.doShape || _ctrl.doFlux) {
257 }
258}
259
261 float normalization = 0.5f * std::erfc(-data / std::sqrt(2.0f * variance));
262 if (!(normalization > 0)) {
263 // avoid division by zero; we know the limit at data << -sigma -> 0.
264 return 0.0;
265 }
266 return data + (std::sqrt(0.5f * variance / boost::math::constants::pi<float>()) *
267 std::exp(-0.5f * (data * data) / variance) / normalization);
268}
269
271 return (std::sqrt(2.0f * variance / boost::math::constants::pi<float>()) *
272 std::exp(-0.5f * (mu * mu) / variance)) -
273 mu * std::erfc(mu / std::sqrt(2.0f * variance));
274}
275
276void BlendednessAlgorithm::_measureMoments(afw::image::MaskedImage<float> const& image,
278 afw::table::Key<double> const& instFluxRawKey,
279 afw::table::Key<double> const& instFluxAbsKey,
280 ShapeResultKey const& _shapeRawKey,
281 ShapeResultKey const& _shapeAbsKey) const {
282 if (_ctrl.doFlux || _ctrl.doShape) {
283 if (!child.getTable()->getCentroidSlot().getMeasKey().isValid()) {
285 "Centroid Key must be defined to measure the blendedness instFlux");
286 }
287 }
288 if (_ctrl.doShape) {
289 if (!child.getTable()->getCentroidSlot().getMeasKey().isValid()) {
291 "Shape Key must be defined to measure the blendedness shape");
292 }
293 }
294 if (_ctrl.doShape || _ctrl.doFlux) {
295 bool fatal = false;
296 if (child.getTable()->getCentroidSlot().getFlagKey().isValid()) {
297 if (child.getCentroidFlag()) {
298 // don't set general flag, because even a failed centroid should
299 // just fall back to the peak, and that should be fine for this
300 // measurement.
301 _flagHandler.setValue(child, NO_CENTROID.number, true);
302 }
303 }
304 if (child.getTable()->getShapeSlot().getFlagKey().isValid()) {
305 if (child.getShapeFlag()) {
306 _flagHandler.setValue(child, NO_SHAPE.number, true);
307 _flagHandler.setValue(child, FAILURE.number, true);
308 }
309 }
310 if (!(child.getShape().getDeterminant() >= 0.0)) {
311 // shape flag should have been set already, but we're paranoid
312 _flagHandler.setValue(child, NO_SHAPE.number, true);
313 _flagHandler.setValue(child, FAILURE.number, true);
314 fatal = true;
315 }
316 if (!(std::isfinite(child.getX()) && std::isfinite(child.getY()))) {
317 // shape flag should have been set already, but we're paranoid
318 _flagHandler.setValue(child, NO_CENTROID.number, true);
319 _flagHandler.setValue(child, FAILURE.number, true);
320 fatal = true;
321 }
322 if (fatal) return;
323 }
324
325 if (_ctrl.doShape) {
326 ShapeAccumulator accumulatorRaw;
327 ShapeAccumulator accumulatorAbs;
328 computeMoments(image, child.getCentroid(), child.getShape(), _ctrl.nSigmaWeightMax, accumulatorRaw,
329 accumulatorAbs);
330 if (_ctrl.doFlux) {
331 child.set(instFluxRawKey, accumulatorRaw.getFlux());
332 child.set(instFluxAbsKey, std::max(accumulatorAbs.getFlux(), 0.0));
333 }
334 _shapeRawKey.set(child, accumulatorRaw.getShape());
335 _shapeAbsKey.set(child, accumulatorAbs.getShape());
336 } else if (_ctrl.doFlux) {
337 FluxAccumulator accumulatorRaw;
338 FluxAccumulator accumulatorAbs;
339 computeMoments(image, child.getCentroid(), child.getShape(), _ctrl.nSigmaWeightMax, accumulatorRaw,
340 accumulatorAbs);
341 child.set(instFluxRawKey, accumulatorRaw.getFlux());
342 child.set(instFluxAbsKey, std::max(accumulatorAbs.getFlux(), 0.0));
343 }
344}
345
347 afw::table::SourceRecord& child) const {
348 _measureMoments(image, child, _instFluxChildRaw, _instFluxChildAbs, _shapeChildRaw, _shapeChildAbs);
349}
350
352 afw::table::SourceRecord& child) const {
353 if (_ctrl.doOld) {
354 child.set(_old, computeOldBlendedness(child.getFootprint(), *image.getImage()));
355 }
356 _measureMoments(image, child, _instFluxParentRaw, _instFluxParentAbs, _shapeParentRaw, _shapeParentAbs);
357 if (_ctrl.doFlux) {
358 child.set(_raw, 1.0 - child.get(_instFluxChildRaw) / child.get(_instFluxParentRaw));
359 child.set(_abs, 1.0 - child.get(_instFluxChildAbs) / child.get(_instFluxParentAbs));
360 if (child.get(_instFluxParentAbs) == 0.0) {
361 // We can get NaNs in the absolute measure if both parent and child have only negative
362 // biased-corrected instFluxes (which we clip to zero). We can't really recover from this,
363 // so we should set the flag.
364 _flagHandler.setValue(child, FAILURE.number, true);
365 }
366 }
367}
368
369} // namespace base
370} // namespace meas
371} // namespace lsst
table::Key< std::string > name
Definition Amplifier.cc:116
AmpInfoBoxKey bbox
Definition Amplifier.cc:117
char * data
Definition BaseRecord.cc:61
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
afw::table::Key< afw::table::Array< VariancePixelT > > variance
int y
Definition SpanSet.cc:48
table::Key< int > transform
SpanPixelIterator Iterator
An iterator over points in the Span.
Definition Span.h:50
std::vector< Span >::const_iterator Iterator
Iterator type used by begin() and end().
Definition PixelRegion.h:50
double getDeterminant() const
Return the determinant of the matrix representation.
Definition Quadrupole.h:83
int getX0() const
Return the image's column-origin.
Definition ImageBase.h:306
typename _const_view_t::x_iterator const_x_iterator
A const iterator for traversing the pixels in a row.
Definition ImageBase.h:141
lsst::geom::Box2I getBBox(ImageOrigin origin=PARENT) const
Definition ImageBase.h:445
int getY0() const
Return the image's row-origin.
Definition ImageBase.h:314
x_iterator x_at(int x, int y) const
Return an x_iterator to the point (x, y) in the image.
Definition ImageBase.h:407
A class to represent a 2-dimensional array of pixels.
Definition Image.h:51
A class to manipulate images, masks, and variance as a single object.
Definition MaskedImage.h:74
const_MaskedImageIterator< typename Image::x_iterator, typename Mask::x_iterator, typename Variance::x_iterator > const_x_iterator
A const_iterator to a row of a MaskedImage.
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition BaseRecord.h:151
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Definition BaseRecord.h:164
A class used as a handle to a particular field in a table.
Definition Key.h:53
Defines the fields and offsets for a table.
Definition Schema.h:51
Record class that contains measurements made on a single exposure.
Definition Source.h:78
bool getShapeFlag() const
Return true if the measurement in the Shape slot failed.
Definition Source.h:590
double getX() const
Return the centroid slot x coordinate.
Definition Source.h:594
double getY() const
Return the centroid slot y coordinate.
Definition Source.h:595
CentroidSlotDefinition::MeasValue getCentroid() const
Get the value of the Centroid slot measurement.
Definition Source.h:570
bool getCentroidFlag() const
Return true if the measurement in the Centroid slot failed.
Definition Source.h:578
std::shared_ptr< SourceTable const > getTable() const
Definition Source.h:102
std::shared_ptr< Footprint > getFootprint() const
Definition Source.h:98
ShapeSlotDefinition::MeasValue getShape() const
Get the value of the Shape slot measurement.
Definition Source.h:582
An integer coordinate rectangle.
Definition Box.h:55
bool contains(Point2I const &point) const noexcept
Return true if the box contains the point.
Definition Box.cc:114
T computeSquaredNorm() const
Return the squared L2 norm of the Extent (x^2 + y^2 + ...).
Definition Extent.h:240
A 2D linear coordinate transformation.
static float computeAbsExpectation(float data, float variance)
Compute the posterior expectation value of the true instrumental flux in a pixel from its (Gaussian) ...
void measureChildPixels(afw::image::MaskedImage< float > const &image, afw::table::SourceRecord &child) const
static FlagDefinitionList const & getFlagDefinitions()
static float computeAbsBias(float mu, float variance)
Compute the bias induced by using the absolute value of a pixel instead of its value.
static FlagDefinition const NO_SHAPE
Definition Blendedness.h:72
static FlagDefinition const NO_CENTROID
Definition Blendedness.h:71
static FlagDefinition const FAILURE
Definition Blendedness.h:70
BlendednessAlgorithm(Control const &ctrl, std::string const &name, afw::table::Schema &schema)
void measureParentPixels(afw::image::MaskedImage< float > const &image, afw::table::SourceRecord &child) const
double nSigmaWeightMax
"Radius factor that sets the maximum extent of the weight function (and hence the " "flux measurement...
Definition Blendedness.h:52
bool doOld
"Whether to compute HeavyFootprint dot products (the old deblend.blendedness parameter)" ;
Definition Blendedness.h:44
bool doShape
"Whether to compute quantities related to the Gaussian-weighted shape" ;
Definition Blendedness.h:48
bool doFlux
"Whether to compute quantities related to the Gaussian-weighted flux" ;
Definition Blendedness.h:46
vector-type utility class to build a collection of FlagDefinitions
Definition FlagHandler.h:60
static FlagHandler addFields(afw::table::Schema &schema, std::string const &prefix, FlagDefinitionList const &flagDefs, FlagDefinitionList const &exclDefs=FlagDefinitionList::getEmptyList())
Add Flag fields to a schema, creating a FlagHandler object to manage them.
void setValue(afw::table::BaseRecord &record, std::size_t i, bool value) const
Set the flag field corresponding to the given flag index.
A FunctorKey for ShapeResult.
static ShapeResultKey addFields(afw::table::Schema &schema, std::string const &name, std::string const &doc, UncertaintyEnum uncertainty, afw::table::CoordinateType coordType=afw::table::CoordinateType::PIXEL)
Add the appropriate fields to a Schema, and return a ShapeResultKey that manages them.
virtual void set(afw::table::BaseRecord &record, ShapeResult const &value) const
Set a ShapeResult in the given record.
Reports errors in the logical structure of the program.
Definition Runtime.h:46
T erfc(T... args)
T exp(T... args)
T isfinite(T... args)
T max(T... args)
T min(T... args)
Point< double, 2 > Point2D
Definition Point.h:324
@ NO_UNCERTAINTY
Algorithm provides no uncertainy information at all.
Definition constants.h:44
UnitVector3d centroid(VertexIterator const begin, VertexIterator const end)
T sqrt(T... args)
double _ww
double _wd
double _w
afw::table::Key< double > weight