LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
ApertureFlux.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 /*
3  * LSST Data Management System
4  * Copyright 2008-2015 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/array.hpp"
27 #include "boost/algorithm/string/replace.hpp"
28 
29 #include "ndarray/eigen.h"
30 
33 #include "lsst/afw/table/Source.h"
36 
37 namespace lsst { namespace meas { namespace base {
38 
39 ApertureFluxControl::ApertureFluxControl() : radii(10), maxSincRadius(10.0), shiftKernel("lanczos5") {
40  // defaults here stolen from HSC pipeline defaults
41  static boost::array<double,10> defaultRadii = {{
42  3.0, 4.5, 6.0, 9.0, 12.0, 17.0, 25.0, 35.0, 50.0, 70.0
43  }};
44  std::copy(defaultRadii.begin(), defaultRadii.end(), radii.begin());
45 }
46 
47 namespace {
48 
49 boost::array<FlagDefinition,ApertureFluxAlgorithm::N_FLAGS> const & getFlagDefinitions() {
50  static boost::array<FlagDefinition,ApertureFluxAlgorithm::N_FLAGS> flagDefs = {{
51  {"flag", "flag set if aperture failed for any reason"},
52  {"flag_apertureTruncated", "flag set if aperture did not fit within the measurement image"},
53  {"flag_sincCoeffsTruncated",
54  "flag set if the full sinc coefficient image for aperture %d did not "
55  "fit within the measurement image"}
56  }};
57  return flagDefs;
58 }
59 
60 } // anonymous
61 
62 std::string ApertureFluxAlgorithm::makeFieldPrefix(std::string const & name, double radius) {
63  std::string prefix = (boost::format("%s_%.1f") % name % radius).str();
64  return boost::replace_all_copy(prefix, ".", "_");
65 }
66 
68  afw::table::Schema & schema, std::string const & prefix, std::string const & doc, bool isSinc
69 ) :
70  fluxKey(FluxResultKey::addFields(schema, prefix, doc)),
71  flags(
72  FlagHandler::addFields(
73  schema, prefix,
74  getFlagDefinitions().begin(),
75  getFlagDefinitions().begin() + (isSinc ? 3 : 2)
76  )
77  )
78 {}
79 
81  Control const & ctrl,
82  std::string const & name,
84  daf::base::PropertySet & metadata
85 
86 ) : _ctrl(ctrl),
87  _centroidExtractor(schema, name)
88 {
89  _keys.reserve(ctrl.radii.size());
90  for (std::size_t i = 0; i < ctrl.radii.size(); ++i) {
91  metadata.add(name + "_radii", ctrl.radii[i]);
92  std::string prefix = ApertureFluxAlgorithm::makeFieldPrefix(name, ctrl.radii[i]);
93  std::string doc = (boost::format("flux within %f-pixel aperture") % ctrl.radii[i]).str();
94  _keys.push_back(Keys(schema, prefix, doc, ctrl.radii[i] <= ctrl.maxSincRadius));
95  }
96 }
97 
99  // This should only get called in the case of completely unexpected failures, so it's not terrible
100  // that we just set the general failure flags for all radii here instead of trying to figure out
101  // which ones we've already done. Any known failure modes are handled inside measure().
102  for (std::size_t i = 0; i < _ctrl.radii.size(); ++i) {
103  _keys[i].flags.handleFailure(measRecord, error);
104  }
105 }
106 
108  Result const & result,
109  afw::table::SourceRecord & record,
110  int index
111 ) const {
112  record.set(_keys[index].fluxKey, result);
113  if (result.getFlag(FAILURE)) {
114  _keys[index].flags.setValue(record, FAILURE, true);
115  }
116  if (result.getFlag(APERTURE_TRUNCATED)) {
117  _keys[index].flags.setValue(record, APERTURE_TRUNCATED, true);
118  }
119  if (result.getFlag(SINC_COEFFS_TRUNCATED)) {
120  _keys[index].flags.setValue(record, SINC_COEFFS_TRUNCATED, true);
121  }
122 }
123 
124 namespace {
125 
126 // Helper function for computeSincFlux get Sinc flux coefficients, and handle cases where the coeff
127 // image needs to be clipped to fit in the measurement image
128 template <typename T>
129 CONST_PTR(afw::image::Image<T>) getSincCoeffs(
130  afw::geom::Box2I const & bbox, // measurement image bbox we need to fit inside
131  afw::geom::ellipses::Ellipse const & ellipse, // ellipse that defines the aperture
132  ApertureFluxAlgorithm::Result & result, // result object where we set flags if we do clip
133  ApertureFluxAlgorithm::Control const & ctrl // configuration
134 ) {
135  CONST_PTR(afw::image::Image<T>) cImage = SincCoeffs<T>::get(ellipse.getCore(), 0.0);
136  cImage = afw::math::offsetImage(
137  *cImage,
138  ellipse.getCenter().getX(),
139  ellipse.getCenter().getY(),
140  ctrl.shiftKernel
141  );
142  if (!bbox.contains(cImage->getBBox())) {
143  // We had to clip out at least part part of the coeff image,
144  // but since that's much larger than the aperture (and close
145  // to zero outside the aperture), it may not be a serious
146  // problem.
148  afw::geom::Box2I overlap = cImage->getBBox();
149  overlap.clip(bbox);
150  if (!overlap.contains(afw::geom::Box2I(ellipse.computeBBox()))) {
151  // The clipping was indeed serious, as we we did have to clip within
152  // the aperture; can't expect any decent answer at this point.
154  result.setFlag(ApertureFluxAlgorithm::FAILURE);
155  }
156  cImage = boost::make_shared< afw::image::Image<T> >(*cImage, overlap);
157  }
158  return cImage;
159 }
160 
161 } // anonymous
162 
163 template <typename T>
165  afw::image::Image<T> const & image,
166  afw::geom::ellipses::Ellipse const & ellipse,
167  Control const & ctrl
168 ) {
169  Result result;
170  CONST_PTR(afw::image::Image<T>) cImage = getSincCoeffs<T>(image.getBBox(), ellipse, result, ctrl);
171  if (result.getFlag(APERTURE_TRUNCATED)) return result;
172  afw::image::Image<T> subImage(image, cImage->getBBox());
173  result.flux = (subImage.getArray().template asEigen<Eigen::ArrayXpr>()
174  * cImage->getArray().template asEigen<Eigen::ArrayXpr>()).sum();
175  return result;
176 }
177 
178 template <typename T>
181  afw::geom::ellipses::Ellipse const & ellipse,
182  Control const & ctrl
183 ) {
184  Result result;
185  CONST_PTR(afw::image::Image<T>) cImage = getSincCoeffs<T>(image.getBBox(), ellipse, result, ctrl);
186  if (result.getFlag(APERTURE_TRUNCATED)) return result;
188  result.flux = (subImage.getImage()->getArray().template asEigen<Eigen::ArrayXpr>()
189  * cImage->getArray().template asEigen<Eigen::ArrayXpr>()).sum();
190  result.fluxSigma = std::sqrt(
191  (subImage.getVariance()->getArray().template asEigen<Eigen::ArrayXpr>().template cast<T>()
192  * cImage->getArray().template asEigen<Eigen::ArrayXpr>().square()).sum()
193  );
194  return result;
195 }
196 
197 template <typename T>
199  afw::image::Image<T> const & image,
200  afw::geom::ellipses::Ellipse const & ellipse,
201  Control const & ctrl
202 ) {
203  Result result;
204  afw::geom::ellipses::PixelRegion region(ellipse); // behaves mostly like a Footprint
205  if (!image.getBBox().contains(region.getBBox())) {
206  result.setFlag(APERTURE_TRUNCATED);
207  result.setFlag(FAILURE);
208  return result;
209  }
210  result.flux = 0;
211  for (
212  afw::geom::ellipses::PixelRegion::Iterator spanIter = region.begin(), spanEnd = region.end();
213  spanIter != spanEnd;
214  ++spanIter
215  ) {
216  typename afw::image::Image<T>::x_iterator pixIter = image.x_at(
217  spanIter->getBeginX() - image.getX0(),
218  spanIter->getY() - image.getY0()
219  );
220  result.flux += std::accumulate(pixIter, pixIter + spanIter->getWidth(), 0.0);
221  }
222  return result;
223 }
224 
225 template <typename T>
228  afw::geom::ellipses::Ellipse const & ellipse,
229  Control const & ctrl
230 ) {
231  Result result;
232  afw::geom::ellipses::PixelRegion region(ellipse); // behaves mostly like a Footprint
233  if (!image.getBBox().contains(region.getBBox())) {
234  result.setFlag(APERTURE_TRUNCATED);
235  result.setFlag(FAILURE);
236  return result;
237  }
238  result.flux = 0.0;
239  result.fluxSigma = 0.0;
240  for (
241  afw::geom::ellipses::PixelRegion::Iterator spanIter = region.begin(), spanEnd = region.end();
242  spanIter != spanEnd;
243  ++spanIter
244  ) {
245  typename afw::image::MaskedImage<T>::Image::x_iterator pixIter = image.getImage()->x_at(
246  spanIter->getBeginX() - image.getX0(),
247  spanIter->getY() - image.getY0()
248  );
249  typename afw::image::MaskedImage<T>::Variance::x_iterator varIter = image.getVariance()->x_at(
250  spanIter->getBeginX() - image.getX0(),
251  spanIter->getY() - image.getY0()
252  );
253  result.flux += std::accumulate(pixIter, pixIter + spanIter->getWidth(), 0.0);
254  // we use this to hold variance as we accumulate...
255  result.fluxSigma += std::accumulate(varIter, varIter + spanIter->getWidth(), 0.0);
256  }
257  result.fluxSigma = std::sqrt(result.fluxSigma); // ...and switch back to sigma here.
258  return result;
259 }
260 
261 template <typename T>
263  afw::image::Image<T> const & image,
264  afw::geom::ellipses::Ellipse const & ellipse,
265  Control const & ctrl
266 ) {
267  return (afw::geom::ellipses::Axes(ellipse.getCore()).getB() <= ctrl.maxSincRadius)
268  ? computeSincFlux(image, ellipse, ctrl)
269  : computeNaiveFlux(image, ellipse, ctrl);
270 }
271 
272 template <typename T>
275  afw::geom::ellipses::Ellipse const & ellipse,
276  Control const & ctrl
277 ) {
278  return (afw::geom::ellipses::Axes(ellipse.getCore()).getB() <= ctrl.maxSincRadius)
279  ? computeSincFlux(image, ellipse, ctrl)
280  : computeNaiveFlux(image, ellipse, ctrl);
281 }
282 #define INSTANTIATE(T) \
283  template \
284  ApertureFluxAlgorithm::Result ApertureFluxAlgorithm::computeFlux( \
285  afw::image::Image<T> const &, \
286  afw::geom::ellipses::Ellipse const &, \
287  Control const & \
288  ); \
289  template \
290  ApertureFluxAlgorithm::Result ApertureFluxAlgorithm::computeFlux( \
291  afw::image::MaskedImage<T> const &, \
292  afw::geom::ellipses::Ellipse const &, \
293  Control const & \
294  ); \
295  template \
296  ApertureFluxAlgorithm::Result ApertureFluxAlgorithm::computeSincFlux( \
297  afw::image::Image<T> const &, \
298  afw::geom::ellipses::Ellipse const &, \
299  Control const & \
300  ); \
301  template \
302  ApertureFluxAlgorithm::Result ApertureFluxAlgorithm::computeSincFlux( \
303  afw::image::MaskedImage<T> const &, \
304  afw::geom::ellipses::Ellipse const &, \
305  Control const & \
306  ); \
307  template \
308  ApertureFluxAlgorithm::Result ApertureFluxAlgorithm::computeNaiveFlux( \
309  afw::image::Image<T> const &, \
310  afw::geom::ellipses::Ellipse const &, \
311  Control const & \
312  ); \
313  template \
314  ApertureFluxAlgorithm::Result ApertureFluxAlgorithm::computeNaiveFlux( \
315  afw::image::MaskedImage<T> const &, \
316  afw::geom::ellipses::Ellipse const &, \
317  Control const & \
318  )
319 
320 INSTANTIATE(float);
321 INSTANTIATE(double);
322 
324  Control const & ctrl,
325  std::string const & name,
326  afw::table::SchemaMapper & mapper
327 ) :
328  BaseTransform(name),
329  _ctrl(ctrl)
330 {
331  for (std::size_t i = 0; i < _ctrl.radii.size(); ++i) {
332  for (auto flag = getFlagDefinitions().begin();
333  flag < getFlagDefinitions().begin() + (_ctrl.radii[i] <= _ctrl.maxSincRadius ? 3 : 2); flag++) {
334  mapper.addMapping(mapper.getInputSchema().find<afw::table::Flag>((boost::format("%s_%s") %
336  flag->name).str()).key);
337  }
340  }
341 }
342 
344  afw::table::SourceCatalog const & inputCatalog,
345  afw::table::BaseCatalog & outputCatalog,
346  afw::image::Wcs const & wcs,
347  afw::image::Calib const & calib
348 ) const {
349  checkCatalogSize(inputCatalog, outputCatalog);
350  std::vector<FluxResultKey> fluxKeys;
351  for (std::size_t i = 0; i < _ctrl.radii.size(); ++i) {
352  fluxKeys.push_back(FluxResultKey(inputCatalog.getSchema()[
354  }
355  afw::table::SourceCatalog::const_iterator inSrc = inputCatalog.begin();
356  afw::table::BaseCatalog::iterator outSrc = outputCatalog.begin();
357  {
358  // While noThrow is in scope, converting a negative flux to a magnitude
359  // returns NaN rather than throwing.
361  for (; inSrc != inputCatalog.end() && outSrc != outputCatalog.end(); ++inSrc, ++outSrc) {
362  for (std::size_t i = 0; i < _ctrl.radii.size(); ++i) {
363  FluxResult fluxResult = fluxKeys[i].get(*inSrc);
364  _magKeys[i].set(*outSrc, calib.getMagnitude(fluxResult.flux, fluxResult.fluxSigma));
365  }
366  }
367  }
368 }
369 
370 }}} // namespace lsst::meas::base
371 
double getMagnitude(double const flux) const
Definition: Calib.cc:391
Defines the fields and offsets for a table.
Definition: Schema.h:46
void setFlag(ApertureFluxAlgorithm::FlagBits bit, bool value=true)
Set the flag value associated with the given bit.
Definition: ApertureFlux.h:254
int getWidth() const
Return the number of columns in the image.
Definition: MaskedImage.h:901
ImagePtr getImage(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s image.
Definition: MaskedImage.h:869
virtual void fail(afw::table::SourceRecord &measRecord, MeasurementError *error=NULL) const
Definition: ApertureFlux.cc:98
void clip(Box2I const &other)
Shrink this to ensure that other.contains(*this).
table::Key< std::string > name
Definition: ApCorrMap.cc:71
Eigen matrix objects that present a view into an ndarray::Array.
Schema & editOutputSchema()
Return a reference to the output schema that allows it to be modified in place.
Definition: SchemaMapper.h:29
A custom container class for records, based on std::vector.
Definition: Catalog.h:94
void checkCatalogSize(afw::table::BaseCatalog const &cat1, afw::table::BaseCatalog const &cat2) const
Ensure that catalogs have the same size.
Definition: Transform.h:101
A mapping between the keys of two Schemas, used to copy data between them.
Definition: SchemaMapper.h:19
#define CONST_PTR(...)
Definition: base.h:47
SelectEigenView< T >::Type copy(Eigen::EigenBase< T > const &other)
Copy an arbitrary Eigen expression into a new EigenView.
Definition: eigen.h:390
virtual void operator()(afw::table::SourceCatalog const &inputCatalog, afw::table::BaseCatalog &outputCatalog, afw::image::Wcs const &wcs, afw::image::Calib const &calib) const
tbl::Key< int > wcs
ApertureFluxAlgorithm(Control const &ctrl, std::string const &name, afw::table::Schema &schema, daf::base::PropertySet &metadata)
Definition: ApertureFlux.cc:80
Implementation of the WCS standard for a any projection.
Definition: Wcs.h:107
Exception to be thrown when a measurement algorithm experiences a known failure mode.
Definition: exceptions.h:48
bool getFlag(ApertureFluxAlgorithm::FlagBits bit) const
Return the flag value associated with the given bit.
Definition: ApertureFlux.h:251
void copyResultToRecord(Result const &result, afw::table::SourceRecord &record, int index) const
boost::enable_if< typename ExpressionTraits< Scalar >::IsScalar, Scalar >::type sum(Scalar const &scalar)
Definition: operators.h:1250
An integer coordinate rectangle.
Definition: Box.h:53
Key< T > addMapping(Key< T > const &inputKey, bool doReplace=false)
Add a new field to the output Schema that is a copy of a field in the input Schema.
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
Keys(afw::table::Schema &schema, std::string const &prefix, std::string const &doc, bool isSinc)
Definition: ApertureFlux.cc:67
VariancePtr getVariance(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s variance.
Definition: MaskedImage.h:890
def error
Definition: log.py:108
int getWidth() const
Return the number of columns in the image.
Definition: Image.h:237
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
Definition: fwd.h:55
Schema getSchema() const
Return the schema associated with the catalog&#39;s table.
Definition: Catalog.h:114
static Result computeFlux(afw::image::Image< T > const &image, afw::geom::ellipses::Ellipse const &ellipse, Control const &ctrl=Control())
geom::Box2I getBBox(ImageOrigin const origin=PARENT) const
Definition: MaskedImage.h:905
An ellipse defined by an arbitrary BaseCore and a center point.
Definition: Ellipse.h:50
double maxSincRadius
&quot;Maximum radius (in pixels) for which the sinc algorithm should be used instead of the &quot; &quot;faster naiv...
Definition: ApertureFlux.h:56
static std::string makeFieldPrefix(std::string const &name, double radius)
Definition: ApertureFlux.cc:62
x_iterator x_at(int x, int y) const
Return an x_iterator to the point (x, y) in the image.
Definition: Image.h:329
Iterator class for CatalogT.
Definition: Catalog.h:34
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:77
static Result computeNaiveFlux(afw::image::Image< T > const &image, afw::geom::ellipses::Ellipse const &ellipse, Control const &ctrl=Control())
SchemaItem< T > find(std::string const &name) const
Find a SchemaItem in the Schema by name.
Flux flux
Measured flux in DN.
Definition: FluxUtilities.h:38
if(width!=gim.getWidth()||height!=gim.getHeight()||x0!=gim.getX0()||y0!=gim.getY0())
Definition: saturated.cc:47
bool contains(Point2I const &point) const
Return true if the box contains the point.
Box2I const & getBBox() const
Definition: PixelRegion.h:45
tbl::Schema schema
A FunctorKey for FluxResult.
Definition: FluxUtilities.h:56
int getY0() const
Definition: Image.h:255
ImageT::Ptr 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:55
static MagResultKey addFields(afw::table::Schema &schema, std::string const &name)
#define INSTANTIATE(T)
An ellipse core for the semimajor/semiminor axis and position angle parametrization (a...
Definition: Axes.h:45
geom::Box2I getBBox(ImageOrigin origin=PARENT) const
Definition: Image.h:377
_view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: Image.h:151
Class for storing generic metadata.
Definition: PropertySet.h:82
static Result computeSincFlux(afw::image::Image< T > const &image, afw::geom::ellipses::Ellipse const &ellipse, Control const &ctrl=Control())
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Definition: BaseRecord.h:136
Record class that contains measurements made on a single exposure.
Definition: Source.h:81
std::vector< MagResultKey > _magKeys
Definition: ApertureFlux.h:288
BaseCore const & getCore() const
Return the ellipse core.
Definition: Ellipse.h:75
std::vector< double > radii
&quot;Radius (in pixels) of apertures.&quot; ;
Definition: ApertureFlux.h:50
void add(std::string const &name, T const &value)
Definition: PropertySet.cc:619
FluxErrElement fluxSigma
1-Sigma error (sqrt of variance) on flux in DN.
Definition: FluxUtilities.h:39
Schema const getInputSchema() const
Return the input schema (copy-on-write).
Definition: SchemaMapper.h:23
A class to represent a 2-dimensional array of pixels.
Definition: PSF.h:43
ApertureFluxTransform(Control const &ctrl, std::string const &name, afw::table::SchemaMapper &mapper)
A reusable result struct for flux measurements.
Definition: FluxUtilities.h:37
int getX0() const
Definition: Image.h:247
SafeCentroidExtractor _centroidExtractor
Definition: ApertureFlux.h:228