LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
KronPhotometry.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 /*
3  * LSST Data Management System
4  * Copyright 2008-2015 LSST Corporation.
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 #include <cmath>
26 #include <functional>
27 #include "boost/algorithm/string.hpp"
28 #include "boost/math/constants/constants.hpp"
29 #include "lsst/pex/exceptions.h"
30 #include "lsst/geom/Point.h"
31 #include "lsst/geom/Box.h"
32 #include "lsst/afw/geom/SpanSet.h"
34 #include "lsst/afw/table/Source.h"
38 #include "lsst/afw/detection/Psf.h"
40 #include "lsst/afw/geom/ellipses.h"
41 #include "lsst/meas/base.h"
43 
45 
46 namespace lsst {
47 namespace meas {
48 namespace extensions {
49 namespace photometryKron {
50 
51 namespace {
52 base::FlagDefinitionList flagDefinitions;
53 } // end anonymous
54 
55 base::FlagDefinition const KronFluxAlgorithm::FAILURE = flagDefinitions.addFailureFlag( "general failure flag, set if anything went wrong");
56 base::FlagDefinition const KronFluxAlgorithm::EDGE = flagDefinitions.add("flag_edge", "bad measurement due to image edge");
57 base::FlagDefinition const KronFluxAlgorithm::BAD_SHAPE_NO_PSF = flagDefinitions.add("flag_bad_shape_no_psf", "bad shape and no PSF");
58 base::FlagDefinition const KronFluxAlgorithm::NO_MINIMUM_RADIUS = flagDefinitions.add("flag_no_minimum_radius", "minimum radius could not enforced: no minimum value or PSF");
59 base::FlagDefinition const KronFluxAlgorithm::NO_FALLBACK_RADIUS = flagDefinitions.add("flag_no_fallback_radius", "no minimum radius and no PSF provided");
60 base::FlagDefinition const KronFluxAlgorithm::BAD_RADIUS = flagDefinitions.add("flag_bad_radius", "bad Kron radius");
61 base::FlagDefinition const KronFluxAlgorithm::USED_MINIMUM_RADIUS = flagDefinitions.add("flag_used_minimum_radius", "used the minimum radius for the Kron aperture");
62 base::FlagDefinition const KronFluxAlgorithm::USED_PSF_RADIUS = flagDefinitions.add("flag_used_psf_radius", "used the PSF Kron radius for the Kron aperture");
63 base::FlagDefinition const KronFluxAlgorithm::SMALL_RADIUS = flagDefinitions.add("flag_small_radius", "measured Kron radius was smaller than that of the PSF");
64 base::FlagDefinition const KronFluxAlgorithm::BAD_SHAPE = flagDefinitions.add("flag_bad_shape", "shape for measuring Kron radius is bad; used PSF shape");
65 
67  return flagDefinitions;
68 }
69 
72 
73 namespace {
74 
75 template <typename MaskedImageT>
76  class FootprintFlux {
77 public:
78  explicit FootprintFlux() : _sum(0.0), _sumVar(0.0) {}
79 
81  void reset() {
82  _sum = _sumVar = 0.0;
83  }
84  void reset(afw::detection::Footprint const&) {}
85 
87  void operator()(geom::Point2I const & pos,
88  typename MaskedImageT::Image::Pixel const & ival,
89  typename MaskedImageT::Variance::Pixel const & vval) {
90  _sum += ival;
91  _sumVar += vval;
92  }
93 
95  double getSum() const { return _sum; }
96 
98  double getSumVar() const { return _sumVar; }
99 
100 private:
101  double _sum;
102  double _sumVar;
103 };
104 
105 /************************************************************************************************************/
116 template <typename MaskedImageT, typename WeightImageT>
117 class FootprintFindMoment {
118 public:
119  FootprintFindMoment(MaskedImageT const& mimage,
120  geom::Point2D const& center, // center of the object
121  double const ab, // axis ratio
122  double const theta // rotation of ellipse +ve from x axis
123  ) : _xcen(center.getX()), _ycen(center.getY()),
124  _ab(ab),
125  _cosTheta(::cos(theta)),
126  _sinTheta(::sin(theta)),
127  _sum(0.0), _sumR(0.0),
128 #if 0
129  _sumVar(0.0), _sumRVar(0.0),
130 #endif
131  _imageX0(mimage.getX0()), _imageY0(mimage.getY0())
132  {}
133 
135  void reset() {}
136  void reset(afw::detection::Footprint const& foot) {
137  _sum = _sumR = 0.0;
138 #if 0
139  _sumVar = _sumRVar = 0.0;
140 #endif
141 
142  MaskedImageT const& mimage = this->getImage();
143  geom::Box2I const& bbox(foot.getBBox());
144  int const x0 = bbox.getMinX(), y0 = bbox.getMinY(), x1 = bbox.getMaxX(), y1 = bbox.getMaxY();
145 
146  if (x0 < _imageX0 || y0 < _imageY0 ||
147  x1 >= _imageX0 + mimage.getWidth() || y1 >= _imageY0 + mimage.getHeight()) {
149  (boost::format("Footprint %d,%d--%d,%d doesn't fit in image %d,%d--%d,%d")
150  % x0 % y0 % x1 % y1
151  % _imageX0 % _imageY0
152  % (_imageX0 + mimage.getWidth() - 1) % (_imageY0 + mimage.getHeight() - 1)
153  ).str());
154  }
155  }
156 
158  void operator()(geom::Point2I const & pos, typename MaskedImageT::Image::Pixel const & ival) {
159  double x = static_cast<double>(pos.getX());
160  double y = static_cast<double>(pos.getY());
161  double const dx = x - _xcen;
162  double const dy = y - _ycen;
163  double const du = dx*_cosTheta + dy*_sinTheta;
164  double const dv = -dx*_sinTheta + dy*_cosTheta;
165 
166  double r = ::hypot(du, dv*_ab); // ellipsoidal radius
167 #if 1
168  if (::hypot(dx, dy) < 0.5) { // within a pixel of the centre
169  /*
170  * We gain significant precision for flattened Gaussians by treating the central pixel specially
171  *
172  * If the object's centered in the pixel (and has constant surface brightness) we have <r> == eR;
173  * if it's at the corner <r> = 2*eR; we interpolate between these exact results linearily in the
174  * displacement. And then add in quadrature which is also a bit dubious
175  *
176  * We could avoid all these issues by estimating <r> using the same trick as we use for
177  * the sinc fluxes; it's not clear that it's worth it.
178  */
179 
180  double const eR = 0.38259771140356325; // <r> for a single square pixel, about the centre
181  r = ::hypot(r, eR*(1 + ::hypot(dx, dy)/geom::ROOT2));
182  }
183 #endif
184 
185  _sum += ival;
186  _sumR += r*ival;
187 #if 0
188  typename MaskedImageT::Variance::Pixel vval = iloc.variance(0, 0);
189  _sumVar += vval;
190  _sumRVar += r*r*vval;
191 #endif
192  }
193 
195  double getIr() const { return _sumR/_sum; }
196 
197 #if 0
199 // double getIrVar() const { return _sumRVar/_sum - getIr()*getIr(); } // Wrong?
200  double getIrVar() const { return _sumRVar/(_sum*_sum) + _sumVar*_sumR*_sumR/::pow(_sum, 4); }
201 #endif
202 
204  bool getGood() const { return _sum > 0 && _sumR > 0; }
205 
206 private:
207  double const _xcen; // center of object
208  double const _ycen; // center of object
209  double const _ab; // axis ratio
210  double const _cosTheta, _sinTheta; // {cos,sin}(angle from x-axis)
211  double _sum; // sum of I
212  double _sumR; // sum of R*I
213 #if 0
214  double _sumVar; // sum of Var(I)
215  double _sumRVar; // sum of R*R*Var(I)
216 #endif
217  int const _imageX0, _imageY0; // origin of image we're measuring
218 
219 };
220 } // end anonymous namespace
221 
223  afw::geom::ellipses::Axes const& shape,
224  geom::LinearTransform const& transformation,
225  double const radius
226  )
227 {
228  afw::geom::ellipses::Axes axes(shape);
229  axes.scale(radius/axes.getDeterminantRadius());
230  return axes.transform(transformation);
231 }
232 
233 template<typename ImageT>
235  ImageT const& image,
237  geom::Point2D const& center,
238  KronFluxControl const& ctrl
239  )
240 {
241  //
242  // We might smooth the image because this is what SExtractor and Pan-STARRS do. But I don't see much gain
243  //
244  double const sigma = ctrl.smoothingSigma; // Gaussian width of smoothing sigma to apply
245  bool const smoothImage = sigma > 0;
246  int kSize = smoothImage ? 2*int(2*sigma) + 1 : 1;
247  afw::math::GaussianFunction1<afw::math::Kernel::Pixel> gaussFunc(smoothImage ? sigma : 100);
248  afw::math::SeparableKernel kernel(kSize, kSize, gaussFunc, gaussFunc);
249  bool const doNormalize = true, doCopyEdge = false;
250  afw::math::ConvolutionControl convCtrl(doNormalize, doCopyEdge);
251  double radius0 = axes.getDeterminantRadius();
253  float radiusForRadius = std::nanf("");
254  for (int i = 0; i < ctrl.nIterForRadius; ++i) {
255  axes.scale(ctrl.nSigmaForRadius);
256  radiusForRadius = axes.getDeterminantRadius(); // radius we used to estimate R_K
257  //
258  // Build an elliptical Footprint of the proper size
259  //
261  afw::geom::ellipses::Ellipse(axes, center)));
262  geom::Box2I bbox = !smoothImage ?
263  foot.getBBox() :
264  kernel.growBBox(foot.getBBox()); // the smallest bbox needed to convolve with Kernel
265  bbox.clip(image.getBBox());
266  ImageT subImage(image, bbox, afw::image::PARENT, smoothImage);
267  if (smoothImage) {
268  afw::math::convolve(subImage, ImageT(image, bbox, afw::image::PARENT, false), kernel, convCtrl);
269  }
270  //
271  // Find the desired first moment of the elliptical radius, which corresponds to the major axis.
272  //
273  FootprintFindMoment<ImageT, afw::detection::Psf::Image> iRFunctor(
274  subImage, center, axes.getA()/axes.getB(), axes.getTheta()
275  );
276 
277  try {
278  foot.getSpans()->applyFunctor(
279  iRFunctor, *(subImage.getImage()));
281  if (i == 0) {
282  LSST_EXCEPT_ADD(e, "Determining Kron aperture");
283  }
284  break; // use the radius we have
285  }
286 
287  if (!iRFunctor.getGood()) {
288  throw LSST_EXCEPT(BadKronException, "Bad integral defining Kron radius");
289  }
290 
291  radius = iRFunctor.getIr()*sqrt(axes.getB()/axes.getA());
292  if (radius <= radius0) {
293  break;
294  }
295  radius0 = radius;
296 
297  axes.scale(radius/axes.getDeterminantRadius()); // set axes to our current estimate of R_K
298  iRFunctor.reset();
299  }
300 
301  return std::make_shared<KronAperture>(center, axes, radiusForRadius);
302 }
303 
304 // Photometer an image with a particular aperture
305 template<typename ImageT>
307  ImageT const& image, // Image to measure
308  afw::geom::ellipses::Ellipse const& aperture, // Aperture in which to measure
309  double const maxSincRadius // largest radius that we use sinc apertures to measure
310  )
311 {
312  afw::geom::ellipses::Axes const& axes = aperture.getCore();
313  if (axes.getB() > maxSincRadius) {
314  FootprintFlux<ImageT> fluxFunctor;
315  auto spans = afw::geom::SpanSet::fromShape(aperture);
316  spans->applyFunctor(
317  fluxFunctor, *(image.getImage()), *(image.getVariance()));
318  return std::make_pair(fluxFunctor.getSum(), ::sqrt(fluxFunctor.getSumVar()));
319  }
320  try {
321  base::ApertureFluxResult fluxResult = base::ApertureFluxAlgorithm::computeSincFlux<float>(image, aperture);
322  return std::make_pair(fluxResult.instFlux, fluxResult.instFluxErr);
323  } catch(pex::exceptions::LengthError &e) {
324  LSST_EXCEPT_ADD(e, (boost::format("Measuring Kron flux for object at (%.3f, %.3f);"
325  " aperture radius %g,%g theta %g")
326  % aperture.getCenter().getX() % aperture.getCenter().getY()
327  % axes.getA() % axes.getB() % geom::radToDeg(axes.getTheta())).str());
328  throw e;
329  }
330 }
331 
332 
334  CONST_PTR(afw::detection::Psf) const& psf, // PSF to measure
335  geom::Point2D const& center, // Centroid of source on parent image
336  double smoothingSigma=0.0 // Gaussian sigma of smoothing applied
337  )
338 {
339  assert(psf);
340  double const radius = psf->computeShape(center).getDeterminantRadius();
341  // For a Gaussian N(0, sigma^2), the Kron radius is sqrt(pi/2)*sigma
342  return ::sqrt(geom::PI/2)*::hypot(radius, std::max(0.0, smoothingSigma));
343 }
344 
345 template<typename ImageT>
347  ImageT const& image,
348  double const nRadiusForFlux,
349  double const maxSincRadius
350  ) const
351 {
352  afw::geom::ellipses::Axes axes(getAxes()); // Copy of ellipse core, so we can scale
353  axes.scale(nRadiusForFlux);
354  afw::geom::ellipses::Ellipse const ellip(axes, getCenter());
355 
356  return photometer(image, ellip, maxSincRadius);
357 }
358 
359 /************************************************************************************************************/
360 
367  KronFluxControl const & ctrl,
368  std::string const & name,
370  daf::base::PropertySet & metadata
371 ) : _name(name),
372  _ctrl(ctrl),
373  _fluxResultKey(
374  meas::base::FluxResultKey::addFields(schema, name, "flux from Kron Flux algorithm")
375  ),
376  _radiusKey(schema.addField<float>(name + "_radius", "Kron radius (sqrt(a*b))")),
377  _radiusForRadiusKey(schema.addField<float>(name + "_radius_for_radius",
378  "radius used to estimate <radius> (sqrt(a*b))")),
379  _psfRadiusKey(schema.addField<float>(name + "_psf_radius", "Radius of PSF")),
380  _centroidExtractor(schema, name, true)
381 {
383  auto metadataName = name + "_nRadiusForflux";
384  boost::to_upper(metadataName);
385  metadata.add(metadataName, ctrl.nRadiusForFlux);
386 }
387 
389  afw::table::SourceRecord & measRecord,
391 ) const {
392  _flagHandler.handleFailure(measRecord, error);
393 }
394 
395 void KronFluxAlgorithm::_applyAperture(
397  afw::image::Exposure<float> const& exposure,
398  KronAperture const& aperture
399  ) const
400 {
401  double const rad = aperture.getAxes().getDeterminantRadius();
403  throw LSST_EXCEPT(
405  BAD_RADIUS.doc,
407  );
408  }
409 
411  try {
412  result = aperture.measureFlux(exposure.getMaskedImage(), _ctrl.nRadiusForFlux, _ctrl.maxSincRadius);
413  } catch (pex::exceptions::LengthError const& e) {
414  // We hit the edge of the image; there's no reasonable fallback or recovery
415  throw LSST_EXCEPT(
417  EDGE.doc,
418  EDGE.number
419  );
421  throw LSST_EXCEPT(
423  EDGE.doc,
424  EDGE.number
425  );
426  }
427 
428  // set the results in the source object
429  meas::base::FluxResult fluxResult;
430  fluxResult.instFlux = result.first;
431  fluxResult.instFluxErr = result.second;
432  source.set(_fluxResultKey, fluxResult);
433  source.set(_radiusKey, aperture.getAxes().getDeterminantRadius());
434  //
435  // REMINDER: In the old code, the psfFactor is calculated using getPsfFactor,
436  // and the values set for _fluxCorrectionKeys. See old meas_algorithms version.
437 }
438 
439 void KronFluxAlgorithm::_applyForced(
440  afw::table::SourceRecord & source,
441  afw::image::Exposure<float> const & exposure,
442  geom::Point2D const & center,
443  afw::table::SourceRecord const & reference,
444  geom::AffineTransform const & refToMeas
445  ) const
446 {
447  float const radius = reference.get(reference.getSchema().find<float>(_ctrl.refRadiusName).key);
448  KronAperture const aperture(reference, refToMeas, radius);
449  _applyAperture(source, exposure, aperture);
450  if (exposure.getPsf()) {
451  source.set(_psfRadiusKey, calculatePsfKronRadius(exposure.getPsf(), center, _ctrl.smoothingSigma));
452  }
453 }
454 
457  afw::image::Exposure<float> const& exposure
458  ) const {
459  geom::Point2D center = _centroidExtractor(source, _flagHandler);
460 
461  // Did we hit a condition that fundamentally prevented measuring the Kron flux?
462  // Such conditions include hitting the edge of the image and bad input shape, but not low signal-to-noise.
463  bool bad = false;
464 
465  afw::image::MaskedImage<float> const& mimage = exposure.getMaskedImage();
466 
467  double R_K_psf = -1;
468  if (exposure.getPsf()) {
469  R_K_psf = calculatePsfKronRadius(exposure.getPsf(), center, _ctrl.smoothingSigma);
470  }
471 
472  //
473  // Get the shape of the desired aperture
474  //
476  if (!source.getShapeFlag()) {
477  axes = source.getShape();
478  } else {
479  bad = true;
480  if (!exposure.getPsf()) {
481  throw LSST_EXCEPT(
485  );
486  }
487  axes = exposure.getPsf()->computeShape();
488  _flagHandler.setValue(source, BAD_SHAPE.number, true);
489  }
490  if (_ctrl.useFootprintRadius) {
491  afw::geom::ellipses::Axes footprintAxes(source.getFootprint()->getShape());
492  // if the Footprint's a disk of radius R we want footRadius == R.
493  // As <r^2> = R^2/2 for a disk, we need to scale up by sqrt(2)
494  footprintAxes.scale(::sqrt(2));
495 
496  double radius0 = axes.getDeterminantRadius();
497  double const footRadius = footprintAxes.getDeterminantRadius();
498 
499  if (footRadius > radius0*_ctrl.nSigmaForRadius) {
500  radius0 = footRadius/_ctrl.nSigmaForRadius; // we'll scale it up by nSigmaForRadius
501  axes.scale(radius0/axes.getDeterminantRadius());
502  }
503  }
504 
505  PTR(KronAperture) aperture;
506  if (_ctrl.fixed) {
507  aperture.reset(new KronAperture(source));
508  } else {
509  try {
510  aperture = KronAperture::determineRadius(mimage, axes, center, _ctrl);
511  } catch (pex::exceptions::OutOfRangeError& e) {
512  // We hit the edge of the image: no reasonable fallback or recovery possible
513  throw LSST_EXCEPT(
515  EDGE.doc,
516  EDGE.number
517  );
518  } catch (BadKronException& e) {
519  // Not setting bad=true because we only failed due to low S/N
520  aperture = _fallbackRadius(source, R_K_psf, e);
521  } catch(pex::exceptions::Exception& e) {
522  bad = true; // There's something fundamental keeping us from measuring the Kron aperture
523  aperture = _fallbackRadius(source, R_K_psf, e);
524  }
525  }
526 
527  /*
528  * Estimate the minimum acceptable Kron radius as the Kron radius of the PSF or the
529  * provided minimum radius
530  */
531 
532  // Enforce constraints on minimum radius
533  double rad = aperture->getAxes().getDeterminantRadius();
534  if (_ctrl.enforceMinimumRadius) {
535  double newRadius = rad;
536  if (_ctrl.minimumRadius > 0.0) {
537  if (rad < _ctrl.minimumRadius) {
538  newRadius = _ctrl.minimumRadius;
539  _flagHandler.setValue(source, USED_MINIMUM_RADIUS.number, true);
540  }
541  } else if (!exposure.getPsf()) {
542  throw LSST_EXCEPT(
546  );
547  } else if (rad < R_K_psf) {
548  newRadius = R_K_psf;
549  _flagHandler.setValue(source, USED_PSF_RADIUS.number, true);
550  }
551  if (newRadius != rad) {
552  aperture->getAxes().scale(newRadius/rad);
553  _flagHandler.setValue(source, SMALL_RADIUS.number, true); // guilty after all
554  }
555  }
556 
557  _applyAperture(source, exposure, *aperture);
558  source.set(_radiusForRadiusKey, aperture->getRadiusForRadius());
559  source.set(_psfRadiusKey, R_K_psf);
560  if (bad) _flagHandler.setValue(source, FAILURE.number, true);
561 }
562 
564  afw::table::SourceRecord & measRecord,
565  afw::image::Exposure<float> const & exposure,
566  afw::table::SourceRecord const & refRecord,
567  afw::geom::SkyWcs const & refWcs
568  ) const {
569  geom::Point2D center = _centroidExtractor(measRecord, _flagHandler);
570  auto xytransform = afw::geom::makeWcsPairTransform(refWcs, *exposure.getWcs());
571  _applyForced(measRecord, exposure, center, refRecord,
572  linearizeTransform(*xytransform, refRecord.getCentroid())
573  );
574 
575 }
576 
577 
578 PTR(KronAperture) KronFluxAlgorithm::_fallbackRadius(afw::table::SourceRecord& source, double const R_K_psf,
579  pex::exceptions::Exception& exc) const
580 {
581  _flagHandler.setValue(source, BAD_RADIUS.number, true);
582  double newRadius;
583  if (_ctrl.minimumRadius > 0) {
584  newRadius = _ctrl.minimumRadius;
585  _flagHandler.setValue(source, USED_MINIMUM_RADIUS.number, true);
586  } else if (R_K_psf > 0) {
587  newRadius = R_K_psf;
588  _flagHandler.setValue(source, USED_PSF_RADIUS.number, true);
589  } else {
590  throw LSST_EXCEPT(
594  );
595  }
596  PTR(KronAperture) aperture(new KronAperture(source));
597  aperture->getAxes().scale(newRadius/aperture->getAxes().getDeterminantRadius());
598  return aperture;
599 }
600 
601 
602 #define INSTANTIATE(TYPE) \
603 template PTR(KronAperture) KronAperture::determineRadius<afw::image::MaskedImage<TYPE> >( \
604  afw::image::MaskedImage<TYPE> const&, \
605  afw::geom::ellipses::Axes, \
606  geom::Point2D const&, \
607  KronFluxControl const& \
608  ); \
609 template std::pair<double, double> KronAperture::measureFlux<afw::image::MaskedImage<TYPE> >( \
610  afw::image::MaskedImage<TYPE> const&, \
611  double const, \
612  double const \
613  ) const;
614 
615 INSTANTIATE(float);
616 
617 }}}} // namespace lsst::meas::extensions::photometryKron
py::object result
Definition: _schema.cc:430
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
double x
#define LSST_EXCEPTION_TYPE(t, b, c)
Macro used to define new types of exceptions without additional data.
Definition: Exception.h:69
#define LSST_EXCEPT_ADD(e, m)
Add the current location and a message to an existing exception before rethrowing it.
Definition: Exception.h:54
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:50
#define INSTANTIATE(TYPE)
int y
Definition: SpanSet.cc:49
table::Schema schema
Definition: python.h:134
#define PTR(...)
Definition: base.h:41
#define CONST_PTR(...)
A shared pointer to a const object.
Definition: base.h:47
Class to describe the properties of a detected object from an image.
Definition: Footprint.h:63
std::shared_ptr< geom::SpanSet > getSpans() const
Return a shared pointer to the SpanSet.
Definition: Footprint.h:115
lsst::geom::Box2I getBBox() const
Return the Footprint's bounding box.
Definition: Footprint.h:208
A polymorphic base class for representing an image's Point Spread Function.
Definition: Psf.h:76
geom::ellipses::Quadrupole computeShape(lsst::geom::Point2D position=makeNullPoint(), image::Color color=image::Color()) const
Compute the ellipse corresponding to the second moments of the Psf.
Definition: Psf.cc:156
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition: SkyWcs.h:117
static std::shared_ptr< geom::SpanSet > fromShape(int r, Stencil s=Stencil::CIRCLE, lsst::geom::Point2I offset=lsst::geom::Point2I())
Factory function for creating SpanSets from a Stencil.
Definition: SpanSet.cc:689
An ellipse core for the semimajor/semiminor axis and position angle parametrization (a,...
Definition: Axes.h:47
double const getTheta() const
Definition: Axes.h:57
double const getA() const
Definition: Axes.h:51
double const getB() const
Definition: Axes.h:54
void scale(double factor)
Scale the size of the ellipse core by the given factor.
Definition: BaseCore.cc:104
double getDeterminantRadius() const
Return the radius defined as the 4th root of the determinant of the quadrupole matrix.
Definition: BaseCore.cc:118
Transformer transform(lsst::geom::LinearTransform const &transform)
Return the transform that maps the ellipse to the unit circle.
Definition: Transformer.h:116
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
MaskedImageT getMaskedImage()
Return the MaskedImage.
Definition: Exposure.h:228
std::shared_ptr< geom::SkyWcs const > getWcs() const
Definition: Exposure.h:232
std::shared_ptr< lsst::afw::detection::Psf const > getPsf() const
Return the Exposure's Psf object.
Definition: Exposure.h:327
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
Parameters to control convolution.
Definition: ConvolveImage.h:50
lsst::geom::Box2I growBBox(lsst::geom::Box2I const &bbox) const
Given a bounding box for pixels one wishes to compute by convolving an image with this kernel,...
Definition: Kernel.cc:177
A kernel described by a pair of functions: func(x, y) = colFunc(x) * rowFunc(y)
Definition: Kernel.h:861
Defines the fields and offsets for a table.
Definition: Schema.h:50
Record class that contains measurements made on a single exposure.
Definition: Source.h:80
CentroidSlotDefinition::MeasValue getCentroid() const
Get the value of the Centroid slot measurement.
Definition: Source.h:572
Class for storing generic metadata.
Definition: PropertySet.h:67
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 affine coordinate transformation consisting of a linear transformation and an offset.
An integer coordinate rectangle.
Definition: Box.h:55
A 2D linear coordinate transformation.
vector-type utility class to build a collection of FlagDefinitions
Definition: FlagHandler.h:60
void handleFailure(afw::table::BaseRecord &record, MeasurementError const *error=nullptr) const
Handle an expected or unexpected Exception thrown by a measurement algorithm.
Definition: FlagHandler.cc:76
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.
Definition: FlagHandler.cc:37
void setValue(afw::table::BaseRecord &record, std::size_t i, bool value) const
Set the flag field corresponding to the given flag index.
Definition: FlagHandler.h:262
Exception to be thrown when a measurement algorithm experiences a known failure mode.
Definition: exceptions.h:48
std::pair< double, double > measureFlux(ImageT const &image, double const nRadiusForFlux, double const maxSincRadius) const
Photometer within the Kron Aperture on an image.
static boost::shared_ptr< KronAperture > determineRadius(ImageT const &image, afw::geom::ellipses::Axes axes, geom::Point2D const &center, KronFluxControl const &ctrl)
Determine the Kron Aperture from an image.
static afw::geom::ellipses::Axes getKronAxes(afw::geom::ellipses::Axes const &shape, geom::LinearTransform const &transformation, double const radius)
Determine Kron axes from a reference image.
static meas::base::FlagDefinition const BAD_RADIUS
static meas::base::FlagDefinition const USED_MINIMUM_RADIUS
static meas::base::FlagDefinition const FAILURE
static meas::base::FlagDefinition const SMALL_RADIUS
virtual void fail(afw::table::SourceRecord &measRecord, meas::base::MeasurementError *error=NULL) const
Handle an exception thrown by the current algorithm by setting flags in the given record.
virtual void measure(afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure) const
Called to measure a single child source in an image.
static meas::base::FlagDefinition const BAD_SHAPE
static meas::base::FlagDefinitionList const & getFlagDefinitions()
virtual void measureForced(afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure, afw::table::SourceRecord const &refRecord, afw::geom::SkyWcs const &refWcs) const
Called to measure a single child source in an image.
static meas::base::FlagDefinition const EDGE
static meas::base::FlagDefinition const NO_MINIMUM_RADIUS
static meas::base::FlagDefinition const NO_FALLBACK_RADIUS
KronFluxAlgorithm(Control const &ctrl, std::string const &name, afw::table::Schema &schema, daf::base::PropertySet &metadata)
A class that knows how to calculate fluxes using the KRON photometry algorithm.
static meas::base::FlagDefinition const BAD_SHAPE_NO_PSF
static meas::base::FlagDefinition const USED_PSF_RADIUS
double smoothingSigma
"Smooth image with N(0, smoothingSigma^2) Gaussian while estimating R_K" ;
std::string refRadiusName
"Name of field specifying reference Kron radius for forced measurement" ;
int nIterForRadius
"Number of times to iterate when setting the Kron radius" ;
bool fixed
"if true, use existing shape and centroid measurements instead of fitting" ;
bool useFootprintRadius
"Use the Footprint size as part of initial estimate of Kron radius" ;
double minimumRadius
"Minimum Kron radius (if == 0.0 use PSF's Kron radius) if enforceMinimumRadius. " "Also functions as ...
double maxSincRadius
"Largest aperture for which to use the slow, accurate, sinc aperture code" ;
double nRadiusForFlux
"Number of Kron radii for Kron flux" ;
bool enforceMinimumRadius
"If true check that the Kron radius exceeds some minimum" ;
double nSigmaForRadius
"Multiplier of rms size for aperture used to initially estimate the Kron radius" ;
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
Reports attempts to exceed implementation-defined length limits for some classes.
Definition: Runtime.h:76
Reports attempts to access elements outside a valid range of indices.
Definition: Runtime.h:89
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104
T hypot(T... args)
T make_pair(T... args)
T max(T... args)
const char * source()
Source function that allows astChannel to source from a Stream.
Definition: Stream.h:224
std::shared_ptr< TransformPoint2ToPoint2 > makeWcsPairTransform(SkyWcs const &src, SkyWcs const &dst)
A Transform obtained by putting two SkyWcs objects "back to back".
Definition: SkyWcs.cc:151
lsst::geom::AffineTransform linearizeTransform(TransformPoint2ToPoint2 const &original, lsst::geom::Point2D const &inPoint)
Approximate a Transform by its local linearization.
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, ConvolutionControl const &convolutionControl=ConvolutionControl())
Convolve an Image or MaskedImage with a Kernel, setting pixels of an existing output image.
lsst::afw::detection::Footprint Footprint
Definition: Source.h:59
constexpr double PI
The ratio of a circle's circumference to diameter.
Definition: Angle.h:39
constexpr double radToDeg(double x) noexcept
Definition: Angle.h:52
constexpr double ROOT2
Definition: Angle.h:46
std::pair< double, double > photometer(ImageT const &image, afw::geom::ellipses::Ellipse const &aperture, double const maxSincRadius)
double calculatePsfKronRadius(boost::shared_ptr< afw::detection::Psf const > const &psf, geom::Point2D const &center, double smoothingSigma=0.0)
float Pixel
Typedefs to be used for pixel values.
Definition: common.h:37
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
double sin(Angle const &a)
Definition: Angle.h:102
double cos(Angle const &a)
Definition: Angle.h:103
A base class for image defects.
T nanf(T... args)
T pow(T... args)
T quiet_NaN(T... args)
A Result struct for running an aperture flux algorithm with a single radius.
Definition: ApertureFlux.h:217
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 > psf
Definition: Exposure.cc:65