LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
PsfFlux.cc
Go to the documentation of this file.
1
2// -*- lsst-c++ -*-
3/*
4 * LSST Data Management System
5 * Copyright 2008-2015 AURA/LSST.
6 *
7 * This product includes software developed by the
8 * LSST Project (http://www.lsst.org/).
9 *
10 * This program is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 3 of the License, or
13 * (at your option) any later version.
14 *
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the LSST License Statement and
21 * the GNU General Public License along with this program. If not,
22 * see <http://www.lsstcorp.org/LegalNotices/>.
23 */
24
25#include <array>
26#include <cmath>
27
28#include "ndarray/eigen.h"
29
32#include "lsst/log/Log.h"
35
36namespace lsst {
37namespace meas {
38namespace base {
39namespace {
40FlagDefinitionList flagDefinitions;
41} // namespace
42
43FlagDefinition const PsfFluxAlgorithm::FAILURE = flagDefinitions.addFailureFlag();
44FlagDefinition const PsfFluxAlgorithm::NO_GOOD_PIXELS =
45 flagDefinitions.add("flag_noGoodPixels", "not enough non-rejected pixels in data to attempt the fit");
46FlagDefinition const PsfFluxAlgorithm::EDGE = flagDefinitions.add(
47 "flag_edge", "object was too close to the edge of the image to use the full PSF model");
48
49FlagDefinitionList const& PsfFluxAlgorithm::getFlagDefinitions() { return flagDefinitions; }
50
51namespace {} // namespace
52
54 std::string const& logName)
55 : _ctrl(ctrl),
56 _instFluxResultKey(FluxResultKey::addFields(
57 schema, name, "instFlux derived from linear least-squares fit of PSF model")),
58 _areaKey(schema.addField<float>(name + "_area", "effective area of PSF", "pixel")),
59 _chi2Key(schema.addField<float>(name + "_chi2", "chi2 of the fitted PSF")),
60 _npixelsKey(schema.addField<int>(name + "_npixels",
61 "number of pixels that were included in the PSF fit", "pixel")),
62 _centroidExtractor(schema, name) {
63 _logName = logName.size() ? logName : name;
64 _flagHandler = FlagHandler::addFields(schema, name, getFlagDefinitions());
65}
66
68 afw::image::Exposure<float> const& exposure) const {
69 std::shared_ptr<afw::detection::Psf const> psf = exposure.getPsf();
70 if (!psf) {
71 LOGL_ERROR(getLogName(), "PsfFlux: no psf attached to exposure");
72 throw LSST_EXCEPT(FatalAlgorithmError, "PsfFlux algorithm requires a Psf with every exposure");
73 }
74 geom::Point2D position = _centroidExtractor(measRecord, _flagHandler);
75 std::shared_ptr<afw::detection::Psf::Image> psfImage = psf->computeImage(position);
76 geom::Box2I fitBBox = psfImage->getBBox();
77 fitBBox.clip(exposure.getBBox());
78 if (fitBBox != psfImage->getBBox()) {
79 _flagHandler.setValue(measRecord, FAILURE.number,
80 true); // if we had a suspect flag, we'd set that instead
81 _flagHandler.setValue(measRecord, EDGE.number, true);
82 }
83 auto fitRegionSpans = std::make_shared<afw::geom::SpanSet>(fitBBox);
84 afw::detection::Footprint fitRegion(fitRegionSpans);
85 if (!_ctrl.badMaskPlanes.empty()) {
86 afw::image::MaskPixel badBits = 0x0;
87 for (std::vector<std::string>::const_iterator i = _ctrl.badMaskPlanes.begin();
88 i != _ctrl.badMaskPlanes.end(); ++i) {
89 badBits |= exposure.getMaskedImage().getMask()->getPlaneBitMask(*i);
90 }
91 fitRegion.setSpans(fitRegion.getSpans()
92 ->intersectNot(*exposure.getMaskedImage().getMask(), badBits)
93 ->clippedTo(exposure.getMaskedImage().getMask()->getBBox()));
94 }
95 if (fitRegion.getArea() == 0) {
96 _flagHandler.setValue(measRecord, NO_GOOD_PIXELS.number, true);
97 _flagHandler.setValue(measRecord, FAILURE.number, true);
98 return;
99 }
100 typedef afw::detection::Psf::Pixel PsfPixel;
101 // SpanSet::flatten returns a new ndarray::Array, which must stay in scope
102 // while we use an Eigen::Map view of it
103 auto modelNdArray = fitRegion.getSpans()->flatten(psfImage->getArray(), psfImage->getXY0());
104 auto dataNdArray = fitRegion.getSpans()->flatten(exposure.getMaskedImage().getImage()->getArray(),
105 exposure.getXY0());
106 auto varianceNdArray = fitRegion.getSpans()->flatten(exposure.getMaskedImage().getVariance()->getArray(),
107 exposure.getXY0());
108 auto model = ndarray::asEigenMatrix(modelNdArray);
109 auto data = ndarray::asEigenMatrix(dataNdArray);
110 auto variance = ndarray::asEigenMatrix(varianceNdArray);
111 PsfPixel alpha = model.squaredNorm();
113 result.instFlux = model.dot(data.cast<PsfPixel>()) / alpha;
114 // If we're not using per-pixel weights to compute the instFlux, we'll still want to compute the
115 // variance as if we had, so we'll apply the weights to the model now, and update alpha.
116 result.instFluxErr = std::sqrt(model.array().square().matrix().dot(variance.cast<PsfPixel>())) / alpha;
117 measRecord.set(_areaKey, model.sum() / alpha);
118 measRecord.set(_npixelsKey, fitRegion.getSpans()->getArea());
119 if (!std::isfinite(result.instFlux) || !std::isfinite(result.instFluxErr)) {
120 throw LSST_EXCEPT(PixelValueError, "Invalid pixel value detected in image.");
121 }
122 measRecord.set(_instFluxResultKey, result);
123 auto chi2 = ((data.cast<PsfPixel>() - result.instFlux * model).array().square() /
124 variance.cast<PsfPixel>().array())
125 .sum();
126 measRecord.set(_chi2Key, chi2);
127}
128
130 _flagHandler.handleFailure(measRecord, error);
131}
132
135 : FluxTransform{name, mapper} {
136 for (std::size_t i = 0; i < PsfFluxAlgorithm::getFlagDefinitions().size(); i++) {
138 if (flag == PsfFluxAlgorithm::FAILURE) continue;
139 if (mapper.getInputSchema().getNames().count(mapper.getInputSchema().join(name, flag.name)) == 0)
140 continue;
142 mapper.getInputSchema().find<afw::table::Flag>(name + "_" + flag.name).key;
143 mapper.addMapping(key);
144 }
145}
146
147} // namespace base
148} // namespace meas
149} // namespace lsst
table::Key< std::string > name
Definition Amplifier.cc:116
char * data
Definition BaseRecord.cc:61
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
afw::table::Key< afw::table::Array< VariancePixelT > > variance
LSST DM logging module built on log4cxx.
#define LOGL_ERROR(logger, message...)
Log a error-level message using a varargs/printf style interface.
Definition Log.h:563
SchemaMapper * mapper
T begin(T... args)
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
std::size_t getArea() const
Return the number of pixels in this Footprint.
Definition Footprint.h:173
void setSpans(std::shared_ptr< geom::SpanSet > otherSpanSet)
Sets the shared pointer to the SpanSet in the Footprint.
Definition Footprint.cc:45
math::Kernel::Pixel Pixel
Pixel type of Image returned by computeImage.
Definition Psf.h:82
A class to contain the data, WCS, and other information needed to describe an image of the sky.
Definition Exposure.h:72
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
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.
Record class that contains measurements made on a single exposure.
Definition Source.h:78
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
std::string getLogName() const
Definition Algorithm.h:66
Exception to be thrown when a measurement algorithm experiences a fatal error.
Definition exceptions.h:76
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
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.
A FunctorKey for FluxResult.
Base for instFlux measurement transformations.
Exception to be thrown when a measurement algorithm experiences a known failure mode.
Definition exceptions.h:48
Exception to be thrown when a measurement algorithm encounters a NaN or infinite pixel.
Definition exceptions.h:83
virtual void measure(afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure) const
Called to measure a single child source in an image.
Definition PsfFlux.cc:67
static FlagDefinition const FAILURE
Definition PsfFlux.h:73
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 PsfFlux.cc:129
static FlagDefinition const NO_GOOD_PIXELS
Definition PsfFlux.h:74
static FlagDefinition const EDGE
Definition PsfFlux.h:75
static FlagDefinitionList const & getFlagDefinitions()
Definition PsfFlux.cc:49
PsfFluxAlgorithm(Control const &ctrl, std::string const &name, afw::table::Schema &schema, std::string const &logName="")
Definition PsfFlux.cc:53
A C++ control class to handle PsfFluxAlgorithm's configuration.
Definition PsfFlux.h:48
std::vector< std::string > badMaskPlanes
"Mask planes that indicate pixels that should be excluded from the fit" ;
Definition PsfFlux.h:51
PsfFluxTransform(Control const &ctrl, std::string const &name, afw::table::SchemaMapper &mapper)
Definition PsfFlux.cc:133
T empty(T... args)
T end(T... args)
T isfinite(T... args)
T size(T... args)
T sqrt(T... args)
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.