LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
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 
32 #include "lsst/afw/table/Source.h"
35 
36 namespace lsst {
37 namespace meas {
38 namespace base {
39 namespace {
40 FlagDefinitionList flagDefinitions;
41 } // namespace
42 
43 FlagDefinition const ApertureFluxAlgorithm::FAILURE = flagDefinitions.addFailureFlag();
44 FlagDefinition const ApertureFluxAlgorithm::APERTURE_TRUNCATED =
45  flagDefinitions.add("flag_apertureTruncated", "aperture did not fit within measurement image");
46 FlagDefinition const ApertureFluxAlgorithm::SINC_COEFFS_TRUNCATED = flagDefinitions.add(
47  "flag_sincCoeffsTruncated", "full sinc coefficient image did not fit within measurement image");
48 
50 
51 ApertureFluxControl::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 
62 ApertureFluxAlgorithm::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 
110 namespace {
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
114 template <typename T>
116 getSincCoeffs(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 ) {
121 std::shared_ptr<afw::image::Image<T> const> cImage = SincCoeffs<T>::get(ellipse.getCore(), 0.0);
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.
136  result.setFlag(ApertureFluxAlgorithm::FAILURE.number);
137  }
138  cImage = std::make_shared<afw::image::Image<T> >(*cImage, overlap);
139  }
140  return cImage;
141 }
142 
143 } // namespace
144 
145 template <typename T>
147  afw::image::Image<T> const &image, afw::geom::ellipses::Ellipse const &ellipse, Control const &ctrl) {
148  Result result;
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 
157 template <typename T>
160  Control const &ctrl) {
161  Result result;
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;
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 
175 template <typename T>
177  afw::image::Image<T> const &image, afw::geom::ellipses::Ellipse const &ellipse, Control const &ctrl) {
178  Result result;
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 
195 template <typename T>
198  Control const &ctrl) {
199  Result result;
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 
222 template <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 
231 template <typename T>
233  afw::geom::ellipses::Ellipse const &ellipse,
234  Control const &ctrl) {
235  return (afw::geom::ellipses::Axes(ellipse.getCore()).getB() <= ctrl.maxSincRadius)
236  ? computeSincFlux(image, ellipse, ctrl)
237  : computeNaiveFlux(image, ellipse, ctrl);
238 }
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 
253 INSTANTIATE(float);
254 INSTANTIATE(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
table::Key< std::string > name
Definition: Amplifier.cc:116
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
std::string prefix
Definition: SchemaMapper.cc:72
SchemaMapper * mapper
Definition: SchemaMapper.cc:71
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:66
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
BaseCore const & getCore() const
Return the ellipse core.
Definition: Ellipse.h:71
lsst::geom::Box2D computeBBox() const
Return the bounding box of the ellipse.
Definition: Ellipse.cc:55
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:73
int getWidth() const
Return the number of columns in the image.
Definition: MaskedImage.h:1054
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage's variance.
Definition: MaskedImage.h:1051
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
Definition: MaskedImage.h:1018
The photometric calibration of an exposure.
Definition: PhotoCalib.h:114
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Definition: BaseRecord.h:164
Iterator class for CatalogT.
Definition: Catalog.h:40
iterator begin()
Iterator access.
Definition: Catalog.h:401
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.
Definition: SchemaMapper.h:21
typename Base::const_iterator const_iterator
Definition: SortedCatalog.h:50
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.
Definition: ApertureFlux.h:80
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.
Definition: ApertureFlux.cc:71
void copyResultToRecord(Result const &result, afw::table::SourceRecord &record, int index) const
Definition: ApertureFlux.cc:96
static FlagDefinitionList const & getFlagDefinitions()
Definition: ApertureFlux.cc:49
static FlagDefinition const FAILURE
Definition: ApertureFlux.h:85
static FlagDefinition const SINC_COEFFS_TRUNCATED
Definition: ApertureFlux.h:87
static std::string makeFieldPrefix(std::string const &name, double radius)
Construct an appropriate prefix for table fields.
Definition: ApertureFlux.cc:57
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
Definition: ApertureFlux.h:86
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.
Definition: ApertureFlux.cc:87
Configuration object for multiple-aperture flux algorithms.
Definition: ApertureFlux.h:49
std::string shiftKernel
"Warping kernel used to shift Sinc photometry coefficients to different center positions" ;
Definition: ApertureFlux.h:62
double maxSincRadius
"Maximum radius (in pixels) for which the sinc algorithm should be used instead of the " "faster naiv...
Definition: ApertureFlux.h:58
std::vector< double > radii
"Radius (in pixels) of apertures." ;
Definition: ApertureFlux.h:53
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
Definition: FlagHandler.h:125
Utility class for handling flag fields that indicate the failure modes of an algorithm.
Definition: FlagHandler.h:148
A FunctorKey for FluxResult.
Definition: FluxUtilities.h:59
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.
Definition: SincCoeffs.cc:505
T copy(T... args)
T end(T... args)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
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
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
T push_back(T... args)
T reserve(T... args)
T size(T... args)
T sqrt(T... args)
#define INSTANTIATE(T)
A Result struct for running an aperture flux algorithm with a single radius.
Definition: ApertureFlux.h:217
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.
Definition: FluxUtilities.h:41
meas::base::Flux instFlux
Measured instFlux in DN.
Definition: FluxUtilities.h:42
meas::base::FluxErrElement instFluxErr
Standard deviation of instFlux in DN.
Definition: FluxUtilities.h:43
Key< int > photoCalib
Definition: Exposure.cc:67