27#include "boost/algorithm/string.hpp"
28#include "boost/math/constants/constants.hpp"
49namespace photometryKron {
52base::FlagDefinitionList flagDefinitions;
55base::FlagDefinition
const KronFluxAlgorithm::FAILURE = flagDefinitions.addFailureFlag(
"general failure flag, set if anything went wrong");
56base::FlagDefinition
const KronFluxAlgorithm::EDGE = flagDefinitions.add(
"flag_edge",
"bad measurement due to image edge");
63base::FlagDefinition
const KronFluxAlgorithm::SMALL_RADIUS = flagDefinitions.add(
"flag_small_radius",
"measured Kron radius was smaller than that of the PSF");
64base::FlagDefinition
const KronFluxAlgorithm::BAD_SHAPE = flagDefinitions.add(
"flag_bad_shape",
"shape for measuring Kron radius is bad; used PSF shape");
67 return flagDefinitions;
75template <
typename MaskedImageT>
78 explicit FootprintFlux() : _sum(0.0), _sumVar(0.0) {}
88 typename MaskedImageT::Image::Pixel
const & ival,
89 typename MaskedImageT::Variance::Pixel
const & vval) {
95 double getSum()
const {
return _sum; }
98 double getSumVar()
const {
return _sumVar; }
116template <
typename MaskedImageT,
typename WeightImageT>
117class FootprintFindMoment {
119 FootprintFindMoment(MaskedImageT
const& mimage,
123 ) : _xcen(center.getX()), _ycen(center.getY()),
125 _cosTheta(::cos(theta)),
126 _sinTheta(::sin(theta)),
127 _sum(0.0), _sumR(0.0),
129 _sumVar(0.0), _sumRVar(0.0),
131 _imageX0(mimage.getX0()), _imageY0(mimage.getY0())
136 void reset(afw::detection::Footprint
const& foot) {
139 _sumVar = _sumRVar = 0.0;
142 MaskedImageT
const& mimage = this->getImage();
144 int const x0 =
bbox.getMinX(), y0 =
bbox.getMinY(), x1 =
bbox.getMaxX(), y1 =
bbox.getMaxY();
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")
151 % _imageX0 % _imageY0
152 % (_imageX0 + mimage.getWidth() - 1) % (_imageY0 + mimage.getHeight() - 1)
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;
166 double r = ::hypot(du, dv*_ab);
168 if (::hypot(dx, dy) < 0.5) {
180 double const eR = 0.38259771140356325;
188 typename MaskedImageT::Variance::Pixel vval = iloc.variance(0, 0);
190 _sumRVar +=
r*
r*vval;
195 double getIr()
const {
return _sumR/_sum; }
200 double getIrVar()
const {
return _sumRVar/(_sum*_sum) + _sumVar*_sumR*_sumR/::pow(_sum, 4); }
204 bool getGood()
const {
return _sum > 0 && _sumR > 0; }
210 double const _cosTheta, _sinTheta;
217 int const _imageX0, _imageY0;
233template<
typename ImageT>
245 bool const smoothImage =
sigma > 0;
246 int kSize = smoothImage ? 2*int(2*
sigma) + 1 : 1;
249 bool const doNormalize =
true, doCopyEdge =
false;
264 kernel.growBBox(foot.
getBBox());
266 ImageT subImage(
image,
bbox, afw::image::PARENT, smoothImage);
273 FootprintFindMoment<ImageT, afw::detection::Psf::Image> iRFunctor(
279 iRFunctor, *(subImage.getImage()));
287 if (!iRFunctor.getGood()) {
291 radius = iRFunctor.getIr()*sqrt(axes.
getB()/axes.
getA());
292 if (radius <= radius0) {
306 return std::make_shared<KronAperture>(center, axes, radiusForRadius);
310template<
typename ImageT>
314 double const maxSincRadius
318 if (axes.
getB() > maxSincRadius) {
319 FootprintFlux<ImageT> fluxFunctor;
322 fluxFunctor, *(
image.getImage()), *(
image.getVariance()));
323 return std::make_pair(fluxFunctor.getSum(), ::sqrt(fluxFunctor.getSumVar()));
329 LSST_EXCEPT_ADD(e, (boost::format(
"Measuring Kron flux for object at (%.3f, %.3f);"
330 " aperture radius %g,%g theta %g")
341 double smoothingSigma=0.0
345 double const radius =
psf->computeShape(center).getDeterminantRadius();
350template<
typename ImageT>
353 double const nRadiusForFlux,
354 double const maxSincRadius
358 axes.
scale(nRadiusForFlux);
379 meas::base::FluxResultKey::addFields(
schema, name,
"flux from Kron Flux algorithm")
381 _radiusKey(
schema.addField<float>(name +
"_radius",
"Kron radius (sqrt(a*b))")),
382 _radiusForRadiusKey(
schema.addField<float>(name +
"_radius_for_radius",
383 "radius used to estimate <radius> (sqrt(a*b))")),
384 _psfRadiusKey(
schema.addField<float>(name +
"_psf_radius",
"Radius of PSF")),
385 _centroidExtractor(
schema, name, true)
388 auto metadataName = name +
"_nRadiusForflux";
389 boost::to_upper(metadataName);
400void KronFluxAlgorithm::_applyAperture(
434 meas::base::FluxResult fluxResult;
435 fluxResult.instFlux =
result.first;
436 fluxResult.instFluxErr =
result.second;
437 source.set(_fluxResultKey, fluxResult);
444void KronFluxAlgorithm::_applyForced(
445 afw::table::SourceRecord & source,
446 afw::image::Exposure<float>
const & exposure,
448 afw::table::SourceRecord
const & reference,
452 float const radius = reference.get(reference.getSchema().find<
float>(_ctrl.
refRadiusName).key);
453 KronAperture
const aperture(reference, refToMeas, radius);
454 _applyAperture(source, exposure, aperture);
455 if (exposure.getPsf()) {
464 geom::Point2D center = _centroidExtractor(source, _flagHandler);
473 if (exposure.getPsf()) {
481 if (!source.getShapeFlag()) {
482 axes = source.getShape();
485 if (!exposure.getPsf()) {
492 axes = exposure.getPsf()->computeShape(exposure.getPsf()->getAveragePosition());
499 footprintAxes.
scale(::sqrt(2));
525 aperture = _fallbackRadius(source, R_K_psf, e);
528 aperture = _fallbackRadius(source, R_K_psf, e);
538 double rad = aperture->getAxes().getDeterminantRadius();
540 double newRadius = rad;
546 }
else if (!exposure.getPsf()) {
552 }
else if (rad < R_K_psf) {
556 if (newRadius != rad) {
557 aperture->getAxes().scale(newRadius/rad);
562 _applyAperture(source, exposure, *aperture);
563 source.set(_radiusForRadiusKey, aperture->getRadiusForRadius());
564 source.set(_psfRadiusKey, R_K_psf);
574 geom::Point2D center = _centroidExtractor(measRecord, _flagHandler);
576 _applyForced(measRecord, exposure, center, refRecord,
577 linearizeTransform(*xytransform, refRecord.
getCentroid())
591 }
else if (R_K_psf > 0) {
607#define INSTANTIATE(TYPE) \
608template std::shared_ptr<KronAperture> KronAperture::determineRadius<afw::image::MaskedImage<TYPE> >( \
609 afw::image::MaskedImage<TYPE> const&, \
610 afw::geom::ellipses::Axes, \
611 geom::Point2D const&, \
612 KronFluxControl const& \
614template std::pair<double, double> KronAperture::measureFlux<afw::image::MaskedImage<TYPE> >( \
615 afw::image::MaskedImage<TYPE> const&, \
#define INSTANTIATE(FROMSYS, TOSYS)
#define LSST_EXCEPTION_TYPE(t, b, c)
Macro used to define new types of exceptions without additional data.
#define LSST_EXCEPT_ADD(e, m)
Add the current location and a message to an existing exception before rethrowing it.
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
afw::table::Key< double > sigma
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
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.
An ellipse core for the semimajor/semiminor axis and position angle parametrization (a,...
double const getTheta() const
double const getA() const
double const getB() const
void scale(double factor)
Scale the size of the ellipse core by the given factor.
double getDeterminantRadius() const
Return the radius defined as the 4th root of the determinant of the quadrupole matrix.
Transformer transform(lsst::geom::LinearTransform const &transform)
Return the transform that maps the ellipse to the unit circle.
An ellipse defined by an arbitrary BaseCore and a center point.
lsst::geom::Point2D const & getCenter() const
Return the center point.
BaseCore const & getCore() const
Return the ellipse core.
A class to contain the data, WCS, and other information needed to describe an image of the sky.
A class to manipulate images, masks, and variance as a single object.
Parameters to control convolution.
A kernel described by a pair of functions: func(x, y) = colFunc(x) * rowFunc(y)
Defines the fields and offsets for a table.
Record class that contains measurements made on a single exposure.
CentroidSlotDefinition::MeasValue getCentroid() const
Get the value of the Centroid slot measurement.
Class for storing generic metadata.
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.
vector-type utility class to build a collection of FlagDefinitions
void handleFailure(afw::table::BaseRecord &record, MeasurementError const *error=nullptr) const
Handle an expected or unexpected Exception thrown by a measurement algorithm.
static FlagHandler addFields(afw::table::Schema &schema, std::string const &prefix, FlagDefinitionList const &flagDefs, FlagDefinitionList const &exclDefs=FlagDefinitionList::getEmptyList())
Add Flag fields to a schema, creating a FlagHandler object to manage them.
void setValue(afw::table::BaseRecord &record, std::size_t i, bool value) const
Set the flag field corresponding to the given flag index.
Exception to be thrown when a measurement algorithm experiences a known failure mode.
std::pair< double, double > measureFlux(ImageT const &image, double const nRadiusForFlux, double const maxSincRadius) const
Photometer within the Kron Aperture on an image.
afw::geom::ellipses::Axes & getAxes()
geom::Point2D const & getCenter() const
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 std::shared_ptr< KronAperture > determineRadius(ImageT const &image, afw::geom::ellipses::Axes axes, geom::Point2D const ¢er, KronFluxControl const &ctrl)
Determine the Kron Aperture from an 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
C++ control object for Kron flux.
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 maxRadius
"Maximum aperture radius in pixels; used to avoid excess memory consumption for faint objects" ;
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.
Reports attempts to exceed implementation-defined length limits for some classes.
Reports attempts to access elements outside a valid range of indices.
Reports errors that are due to events beyond the control of the program.
const char * source()
Source function that allows astChannel to source from a Stream.
std::shared_ptr< TransformPoint2ToPoint2 > makeWcsPairTransform(SkyWcs const &src, SkyWcs const &dst)
A Transform obtained by putting two SkyWcs objects "back to back".
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.
constexpr double radToDeg(double x) noexcept
double constexpr PI
The ratio of a circle's circumference to diameter.
std::pair< double, double > photometer(ImageT const &image, afw::geom::ellipses::Ellipse const &aperture, double const maxSincRadius)
double calculatePsfKronRadius(std::shared_ptr< afw::detection::Psf const > const &psf, geom::Point2D const ¢er, double smoothingSigma=0.0)
A Result struct for running an aperture flux algorithm with a single radius.
meas::base::Flux instFlux
Measured instFlux in DN.
meas::base::FluxErrElement instFluxErr
Standard deviation of instFlux in DN.