LSST Applications g063fba187b+66a50001ff,g0f08755f38+1a22dc2551,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g217e2c1bcf+12e87a5bd8,g246886dfd9+466c7b6c06,g28da252d5a+858b171e04,g2bbee38e9b+404b60ec9b,g2bc492864f+404b60ec9b,g3156d2b45e+6e55a43351,g347aa1857d+404b60ec9b,g35bb328faa+a8ce1bb630,g3a166c0a6a+404b60ec9b,g3e281a1b8c+c5dd892a6c,g414038480c+6b9177ef31,g41af890bb2+70bea58702,g599934f4f4+b8c5400ca5,g781aacb6e4+a8ce1bb630,g7af13505b9+b5b9cefdb8,g80478fca09+c2997882f3,g82479be7b0+8974e6af0f,g858d7b2824+1a22dc2551,g89c8672015+f4add4ffd5,g8f1c07a47a+de51c9b0a5,g9125e01d80+a8ce1bb630,ga5288a1d22+b66f8cf76b,gb58c049af0+d64f4d3760,gc28159a63d+404b60ec9b,gcab2d0539d+66cf1de5d4,gcf0d15dbbd+12cb7e2563,gda6a2b7d83+12cb7e2563,gdaeeff99f8+1711a396fd,ge79ae78c31+404b60ec9b,gef2f8181fd+414189b318,gf0baf85859+c1f95f4921,gf0c06eb49c+1a22dc2551,gfa517265be+1a22dc2551,gfa999e8aa5+17cd334064,v28.0.0.rc2
LSST Data Management Base Package
Loading...
Searching...
No Matches
ApertureFlux.cc
Go to the documentation of this file.
1// -*- lsst-c++ -*-
2/*
3 * LSST Data Management System
4 * Copyright 2008-2016 AURA/LSST.
5 *
6 * This product includes software developed by the
7 * LSST Project (http://www.lsst.org/).
8 *
9 * This program is free software: you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation, either version 3 of the License, or
12 * (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the LSST License Statement and
20 * the GNU General Public License along with this program. If not,
21 * see <http://www.lsstcorp.org/LegalNotices/>.
22 */
23
24#include <numeric>
25
26#include "boost/algorithm/string.hpp"
27
28#include "ndarray/eigen.h"
29
35
36namespace lsst {
37namespace meas {
38namespace base {
39namespace {
40FlagDefinitionList flagDefinitions;
41} // namespace
42
43FlagDefinition const ApertureFluxAlgorithm::FAILURE = flagDefinitions.addFailureFlag();
45 flagDefinitions.add("flag_apertureTruncated", "aperture did not fit within measurement image");
46FlagDefinition const ApertureFluxAlgorithm::SINC_COEFFS_TRUNCATED = flagDefinitions.add(
47 "flag_sincCoeffsTruncated", "full sinc coefficient image did not fit within measurement image");
48
50
51ApertureFluxControl::ApertureFluxControl() : radii(10), maxSincRadius(10.0), shiftKernel("lanczos5") {
52 // defaults here stolen from HSC pipeline defaults
53 static std::array<double, 10> defaultRadii = {{3.0, 4.5, 6.0, 9.0, 12.0, 17.0, 25.0, 35.0, 50.0, 70.0}};
54 std::copy(defaultRadii.begin(), defaultRadii.end(), radii.begin());
55}
56
58 std::string prefix = (boost::format("%s_%.1f") % name % radius).str();
59 return boost::replace_all_copy(prefix, ".", "_");
60}
61
62ApertureFluxAlgorithm::Keys::Keys(afw::table::Schema &schema, std::string const &prefix,
63 std::string const &doc, bool isSinc)
64 : instFluxKey(FluxResultKey::addFields(schema, prefix, doc)),
65 flags(
66 // The exclusion List can either be empty, or constain the sinc coeffs flag
67 FlagHandler::addFields(
68 schema, prefix, ApertureFluxAlgorithm::getFlagDefinitions(),
70
73
74 )
75 : _ctrl(ctrl), _centroidExtractor(schema, name) {
76 _keys.reserve(ctrl.radii.size());
77 for (std::size_t i = 0; i < ctrl.radii.size(); ++i) {
78 std::string upperName(name);
79 boost::to_upper(upperName);
80 metadata.add(upperName + "_RADII", ctrl.radii[i]);
82 std::string doc = (boost::format("instFlux within %f-pixel aperture") % ctrl.radii[i]).str();
83 _keys.push_back(Keys(schema, prefix, doc, ctrl.radii[i] <= ctrl.maxSincRadius));
84 }
85}
86
88 // This should only get called in the case of completely unexpected failures, so it's not terrible
89 // that we just set the general failure flags for all radii here instead of trying to figure out
90 // which ones we've already done. Any known failure modes are handled inside measure().
91 for (std::size_t i = 0; i < _ctrl.radii.size(); ++i) {
92 _keys[i].flags.handleFailure(measRecord, error);
93 }
94}
95
97 int index) const {
98 record.set(_keys[index].instFluxKey, result);
99 if (result.getFlag(FAILURE.number)) {
100 _keys[index].flags.setValue(record, FAILURE.number, true);
101 }
102 if (result.getFlag(APERTURE_TRUNCATED.number)) {
103 _keys[index].flags.setValue(record, APERTURE_TRUNCATED.number, true);
104 }
105 if (result.getFlag(SINC_COEFFS_TRUNCATED.number)) {
106 _keys[index].flags.setValue(record, SINC_COEFFS_TRUNCATED.number, true);
107 }
108}
109
110namespace {
111
112// Helper function for computeSincFlux get Sinc instFlux coefficients, and handle cases where the coeff
113// image needs to be clipped to fit in the measurement image
114template <typename T>
116getSincCoeffs(geom::Box2I const &bbox, // measurement image bbox we need to fit inside
117 afw::geom::ellipses::Ellipse const &ellipse, // ellipse that defines the aperture
118 ApertureFluxAlgorithm::Result &result, // result object where we set flags if we do clip
119 ApertureFluxAlgorithm::Control const &ctrl // configuration
120) {
122 cImage = afw::math::offsetImage(*cImage, ellipse.getCenter().getX(), ellipse.getCenter().getY(),
123 ctrl.shiftKernel);
124 if (!bbox.contains(cImage->getBBox())) {
125 // We had to clip out at least part part of the coeff image,
126 // but since that's much larger than the aperture (and close
127 // to zero outside the aperture), it may not be a serious
128 // problem.
130 geom::Box2I overlap = cImage->getBBox();
131 overlap.clip(bbox);
132 if (!overlap.contains(geom::Box2I(ellipse.computeBBox()))) {
133 // The clipping was indeed serious, as we we did have to clip within
134 // the aperture; can't expect any decent answer at this point.
137 }
138 cImage = std::make_shared<afw::image::Image<T> >(*cImage, overlap);
139 }
140 return cImage;
141}
142
143} // namespace
144
145template <typename T>
147 afw::image::Image<T> const &image, afw::geom::ellipses::Ellipse const &ellipse, Control const &ctrl) {
149 std::shared_ptr<afw::image::Image<T> const> cImage = getSincCoeffs<T>(image.getBBox(), ellipse, result, ctrl);
150 if (result.getFlag(APERTURE_TRUNCATED.number)) return result;
151 afw::image::Image<T> subImage(image, cImage->getBBox());
152 result.instFlux =
153 (ndarray::asEigenArray(subImage.getArray()) * ndarray::asEigenArray(cImage->getArray())).sum();
154 return result;
155}
156
157template <typename T>
160 Control const &ctrl) {
162 std::shared_ptr<afw::image::Image<T> const> cImage = getSincCoeffs<T>(image.getBBox(), ellipse, result, ctrl);
163 if (result.getFlag(APERTURE_TRUNCATED.number)) return result;
164 afw::image::MaskedImage<T> subImage(image, cImage->getBBox(afw::image::PARENT), afw::image::PARENT);
165 result.instFlux = (ndarray::asEigenArray(subImage.getImage()->getArray()) *
166 ndarray::asEigenArray(cImage->getArray()))
167 .sum();
168 result.instFluxErr =
169 std::sqrt((ndarray::asEigenArray(subImage.getVariance()->getArray()).template cast<T>() *
170 ndarray::asEigenArray(cImage->getArray()).square())
171 .sum());
172 return result;
173}
174
175template <typename T>
177 afw::image::Image<T> const &image, afw::geom::ellipses::Ellipse const &ellipse, Control const &ctrl) {
179 afw::geom::ellipses::PixelRegion region(ellipse); // behaves mostly like a Footprint
180 if (!image.getBBox().contains(region.getBBox())) {
182 result.setFlag(FAILURE.number);
183 return result;
184 }
185 result.instFlux = 0;
186 for (afw::geom::ellipses::PixelRegion::Iterator spanIter = region.begin(), spanEnd = region.end();
187 spanIter != spanEnd; ++spanIter) {
188 typename afw::image::Image<T>::x_iterator pixIter =
189 image.x_at(spanIter->getBeginX() - image.getX0(), spanIter->getY() - image.getY0());
190 result.instFlux += std::accumulate(pixIter, pixIter + spanIter->getWidth(), 0.0);
191 }
192 return result;
193}
194
195template <typename T>
198 Control const &ctrl) {
200 afw::geom::ellipses::PixelRegion region(ellipse); // behaves mostly like a Footprint
201 if (!image.getBBox().contains(region.getBBox())) {
203 result.setFlag(FAILURE.number);
204 return result;
205 }
206 result.instFlux = 0.0;
207 result.instFluxErr = 0.0;
208 for (afw::geom::ellipses::PixelRegion::Iterator spanIter = region.begin(), spanEnd = region.end();
209 spanIter != spanEnd; ++spanIter) {
210 typename afw::image::MaskedImage<T>::Image::x_iterator pixIter = image.getImage()->x_at(
211 spanIter->getBeginX() - image.getX0(), spanIter->getY() - image.getY0());
212 typename afw::image::MaskedImage<T>::Variance::x_iterator varIter = image.getVariance()->x_at(
213 spanIter->getBeginX() - image.getX0(), spanIter->getY() - image.getY0());
214 result.instFlux += std::accumulate(pixIter, pixIter + spanIter->getWidth(), 0.0);
215 // we use this to hold variance as we accumulate...
216 result.instFluxErr += std::accumulate(varIter, varIter + spanIter->getWidth(), 0.0);
217 }
218 result.instFluxErr = std::sqrt(result.instFluxErr); // ...and switch back to sigma here.
219 return result;
220}
221
222template <typename T>
224 afw::geom::ellipses::Ellipse const &ellipse,
225 Control const &ctrl) {
226 return (afw::geom::ellipses::Axes(ellipse.getCore()).getB() <= ctrl.maxSincRadius)
227 ? computeSincFlux(image, ellipse, ctrl)
228 : computeNaiveFlux(image, ellipse, ctrl);
229}
230
231template <typename T>
239#define INSTANTIATE(T) \
240 template ApertureFluxAlgorithm::Result ApertureFluxAlgorithm::computeFlux( \
241 afw::image::Image<T> const &, afw::geom::ellipses::Ellipse const &, Control const &); \
242 template ApertureFluxAlgorithm::Result ApertureFluxAlgorithm::computeFlux( \
243 afw::image::MaskedImage<T> const &, afw::geom::ellipses::Ellipse const &, Control const &); \
244 template ApertureFluxAlgorithm::Result ApertureFluxAlgorithm::computeSincFlux( \
245 afw::image::Image<T> const &, afw::geom::ellipses::Ellipse const &, Control const &); \
246 template ApertureFluxAlgorithm::Result ApertureFluxAlgorithm::computeSincFlux( \
247 afw::image::MaskedImage<T> const &, afw::geom::ellipses::Ellipse const &, Control const &); \
248 template ApertureFluxAlgorithm::Result ApertureFluxAlgorithm::computeNaiveFlux( \
249 afw::image::Image<T> const &, afw::geom::ellipses::Ellipse const &, Control const &); \
250 template ApertureFluxAlgorithm::Result ApertureFluxAlgorithm::computeNaiveFlux( \
251 afw::image::MaskedImage<T> const &, afw::geom::ellipses::Ellipse const &, Control const &)
252
253INSTANTIATE(float);
254INSTANTIATE(double);
255
258 : BaseTransform(name), _ctrl(ctrl) {
259 for (std::size_t i = 0; i < _ctrl.radii.size(); ++i) {
262 if (_ctrl.radii[i] > _ctrl.maxSincRadius &&
264 continue;
265 }
266 mapper.addMapping(
267 mapper.getInputSchema()
268 .find<afw::table::Flag>(
269 (boost::format("%s_%s") %
271 .str())
272 .key);
273 }
274 _magKeys.push_back(MagResultKey::addFields(
275 mapper.editOutputSchema(), ApertureFluxAlgorithm::makeFieldPrefix(name, _ctrl.radii[i])));
276 }
277}
278
280 afw::table::BaseCatalog &outputCatalog, afw::geom::SkyWcs const &wcs,
281 afw::image::PhotoCalib const &photoCalib) const {
282 checkCatalogSize(inputCatalog, outputCatalog);
283 std::vector<FluxResultKey> instFluxKeys;
284 for (std::size_t i = 0; i < _ctrl.radii.size(); ++i) {
285 instFluxKeys.push_back(FluxResultKey(
286 inputCatalog.getSchema()[ApertureFluxAlgorithm::makeFieldPrefix(_name, _ctrl.radii[i])]));
287 }
288 afw::table::SourceCatalog::const_iterator inSrc = inputCatalog.begin();
289 afw::table::BaseCatalog::iterator outSrc = outputCatalog.begin();
290 {
291 for (; inSrc != inputCatalog.end() && outSrc != outputCatalog.end(); ++inSrc, ++outSrc) {
292 for (std::size_t i = 0; i < _ctrl.radii.size(); ++i) {
293 FluxResult instFluxResult = instFluxKeys[i].get(*inSrc);
294 _magKeys[i].set(*outSrc, photoCalib.instFluxToMagnitude(instFluxResult.instFlux,
295 instFluxResult.instFluxErr));
296 }
297 }
298 }
299}
300
301} // namespace base
302} // namespace meas
303} // namespace lsst
py::object result
Definition _schema.cc:429
AmpInfoBoxKey bbox
Definition Amplifier.cc:117
#define INSTANTIATE(FROMSYS, TOSYS)
Definition Detector.cc:509
std::string prefix
SchemaMapper * mapper
T accumulate(T... args)
table::Schema schema
Definition python.h:134
T begin(T... args)
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition SkyWcs.h:117
An ellipse core for the semimajor/semiminor axis and position angle parametrization (a,...
Definition Axes.h:47
double const getB() const
Definition Axes.h:54
An ellipse defined by an arbitrary BaseCore and a center point.
Definition Ellipse.h:51
lsst::geom::Point2D const & getCenter() const
Return the center point.
Definition Ellipse.h:62
lsst::geom::Box2D computeBBox() const
Return the bounding box of the ellipse.
Definition Ellipse.cc:55
BaseCore const & getCore() const
Return the ellipse core.
Definition Ellipse.h:71
A pixelized region containing all pixels whose centers are within an Ellipse.
Definition PixelRegion.h:46
lsst::geom::Box2I const & getBBox() const
Return the bounding box of the pixel region.
Definition PixelRegion.h:80
Iterator begin() const
Iterator range over Spans whose pixels are within the Ellipse.
Definition PixelRegion.h:68
std::vector< Span >::const_iterator Iterator
Iterator type used by begin() and end().
Definition PixelRegion.h:50
int getWidth() const
Return the number of columns in the image.
Definition ImageBase.h:294
typename _view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition ImageBase.h:133
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
int getWidth() const
Return the number of columns in the image.
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage's variance.
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
The photometric calibration of an exposure.
Definition PhotoCalib.h:114
Tag types used to declare specialized field types.
Definition misc.h:31
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Definition BaseRecord.h:164
iterator begin()
Iterator access.
Definition Catalog.h:400
Defines the fields and offsets for a table.
Definition Schema.h:51
A mapping between the keys of two Schemas, used to copy data between them.
typename Base::const_iterator const_iterator
Record class that contains measurements made on a single exposure.
Definition Source.h:78
Class for storing generic metadata.
Definition PropertySet.h:66
void add(std::string const &name, T const &value)
Append a single value to the vector of values for a property name (possibly hierarchical).
An integer coordinate rectangle.
Definition Box.h:55
void clip(Box2I const &other) noexcept
Shrink this to ensure that other.contains(*this).
Definition Box.cc:189
bool contains(Point2I const &point) const noexcept
Return true if the box contains the point.
Definition Box.cc:114
Base class for multiple-aperture photometry algorithms.
ApertureFluxAlgorithm(Control const &ctrl, std::string const &name, afw::table::Schema &schema, daf::base::PropertySet &metadata)
Construct the algorithm and add its fields to the given Schema.
void copyResultToRecord(Result const &result, afw::table::SourceRecord &record, int index) const
static FlagDefinitionList const & getFlagDefinitions()
static FlagDefinition const FAILURE
static FlagDefinition const SINC_COEFFS_TRUNCATED
static std::string makeFieldPrefix(std::string const &name, double radius)
Construct an appropriate prefix for table fields.
static Result computeSincFlux(afw::image::Image< T > const &image, afw::geom::ellipses::Ellipse const &ellipse, Control const &ctrl=Control())
Compute the instFlux (and optionally, uncertanties) within an aperture using Sinc photometry.
static Result computeFlux(afw::image::Image< T > const &image, afw::geom::ellipses::Ellipse const &ellipse, Control const &ctrl=Control())
Compute the instFlux (and optionally, uncertanties) within an aperture using the algorithm determined...
static FlagDefinition const APERTURE_TRUNCATED
static Result computeNaiveFlux(afw::image::Image< T > const &image, afw::geom::ellipses::Ellipse const &ellipse, Control const &ctrl=Control())
Compute the instFlux (and optionally, uncertanties) within an aperture using naive photometry.
virtual void fail(afw::table::SourceRecord &measRecord, MeasurementError *error=nullptr) const
Handle an exception thrown by the current algorithm by setting flags in the given record.
Configuration object for multiple-aperture flux algorithms.
std::string shiftKernel
"Warping kernel used to shift Sinc photometry coefficients to different center positions" ;
double maxSincRadius
"Maximum radius (in pixels) for which the sinc algorithm should be used instead of the " "faster naiv...
std::vector< double > radii
"Radius (in pixels) of apertures." ;
ApertureFluxTransform(Control const &ctrl, std::string const &name, afw::table::SchemaMapper &mapper)
virtual void operator()(afw::table::SourceCatalog const &inputCatalog, afw::table::BaseCatalog &outputCatalog, afw::geom::SkyWcs const &wcs, afw::image::PhotoCalib const &photoCalib) const
Abstract base class for all C++ measurement transformations.
Definition Transform.h:86
void checkCatalogSize(afw::table::BaseCatalog const &cat1, afw::table::BaseCatalog const &cat2) const
Ensure that catalogs have the same size.
Definition Transform.h:102
vector-type utility class to build a collection of FlagDefinitions
Definition FlagHandler.h:60
std::size_t size() const
return the current size (number of defined elements) of the collection
Utility class for handling flag fields that indicate the failure modes of an algorithm.
A FunctorKey for FluxResult.
static MagResultKey addFields(afw::table::Schema &schema, std::string const &name)
Add a pair of _mag, _magErr fields to a Schema, and return a MagResultKey that points to them.
Exception to be thrown when a measurement algorithm experiences a known failure mode.
Definition exceptions.h:48
static std::shared_ptr< CoeffT const > get(afw::geom::ellipses::Axes const &outerEllipse, float const innerRadiusFactor=0.0)
Get the coefficients for an aperture.
T copy(T... args)
T end(T... args)
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 push_back(T... args)
T reserve(T... args)
T size(T... args)
T sqrt(T... args)
A Result struct for running an aperture flux algorithm with a single radius.
Simple class used to define and document flags The name and doc constitute the identity of the FlagDe...
Definition FlagHandler.h:40
A reusable result struct for instFlux measurements.
meas::base::Flux instFlux
Measured instFlux in DN.
meas::base::FluxErrElement instFluxErr
Standard deviation of instFlux in DN.