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
PeakLikelihoodFlux.cc
Go to the documentation of this file.
1// -*- lsst-c++ -*-
2/*
3 * LSST Data Management System
4 * Copyright 2008-2013 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 "ndarray/eigen.h"
25
27#include "lsst/geom/Box.h"
28#include "lsst/pex/exceptions.h"
29#include "lsst/afw/geom.h"
30#include "lsst/afw/image.h"
31#include "lsst/afw/math.h"
32
35
36namespace lsst {
37namespace meas {
38namespace base {
39namespace {
40FlagDefinitionList flagDefinitions;
41} // namespace
42
43FlagDefinition const PeakLikelihoodFluxAlgorithm::FAILURE = flagDefinitions.addFailureFlag();
44
46
47namespace {
48
49/************************************************************************************************************/
66class PsfAttributes {
67public:
68 enum Method {
69 ADAPTIVE_MOMENT,
70 FIRST_MOMENT,
71 SECOND_MOMENT,
72 NOISE_EQUIVALENT,
73 BICKERTON
74 };
75
76 PsfAttributes(std::shared_ptr<afw::detection::Psf const> psf, int const iX, int const iY);
77 PsfAttributes(std::shared_ptr<afw::detection::Psf const> psf, geom::Point2I const &cen);
78
79 double computeGaussianWidth(Method how = ADAPTIVE_MOMENT) const;
80 double computeEffectiveArea() const;
81
82private:
84};
85
89PsfAttributes::PsfAttributes(std::shared_ptr<afw::detection::Psf const> psf,
90 int const iX,
91 int const iY
92 ) {
93 // N.b. (iX, iY) are ints so that we know this image is centered in the central pixel of _psfImage
94 _psfImage = psf->computeImage(geom::PointD(iX, iY));
95}
96
100PsfAttributes::PsfAttributes(
102 geom::Point2I const &cen
103 )
104 : // N.b. cen is a PointI so that we know this image is centered in the central pixel of _psfImage
105 _psfImage(psf->computeImage(geom::PointD(cen))) {}
106
111double PsfAttributes::computeEffectiveArea() const {
112 double sum = 0.0;
113 double sumsqr = 0.0;
114 for (int iY = 0; iY != _psfImage->getHeight(); ++iY) {
115 afw::image::Image<double>::x_iterator end = _psfImage->row_end(iY);
116 for (afw::image::Image<double>::x_iterator ptr = _psfImage->row_begin(iY); ptr != end; ++ptr) {
117 sum += *ptr;
118 sumsqr += (*ptr) * (*ptr);
119 }
120 }
121 return sum * sum / sumsqr;
122}
123
124} // end anonymous namespace
125
133template <typename T>
135 afw::image::MaskedImage<T> const &maskedImage,
137 geom::Point2D const &fracShift,
138 geom::Point2I const &parentInd
139 ) {
140 typedef typename afw::image::Exposure<T>::MaskedImageT MaskedImageT;
141 typedef typename afw::image::Image<double> KernelImageT;
142
144
145 if ((std::abs(fracShift[0]) >= 1) || (std::abs(fracShift[1]) >= 1)) {
147 os << "fracShift = " << fracShift << " too large; abs value must be < 1 in both axes";
149 }
150
151 // warping kernels have even dimension and want the peak to the right of center
152 if (fracShift[0] < 0) {
153 warpingKernelPtr->setCtr(warpingKernelPtr->getCtr() + geom::Extent2I(1, 0));
154 }
155 if (fracShift[1] < 0) {
156 warpingKernelPtr->setCtr(warpingKernelPtr->getCtr() + geom::Extent2I(0, 1));
157 }
158 geom::Box2I warpingOverlapBBox(parentInd - geom::Extent2I(warpingKernelPtr->getCtr()),
159 warpingKernelPtr->getDimensions());
160 if (!maskedImage.getBBox().contains(warpingOverlapBBox)) {
162 os << "Warping kernel extends off the edge"
163 << "; kernel bbox = " << warpingOverlapBBox << "; exposure bbox = " << maskedImage.getBBox();
165 }
166 warpingKernelPtr->setKernelParameters(std::make_pair(fracShift[0], fracShift[1]));
167 KernelImageT warpingKernelImage(warpingKernelPtr->getDimensions());
168 warpingKernelPtr->computeImage(warpingKernelImage, true);
169 typename KernelImageT::const_xy_locator const warpingKernelLoc = warpingKernelImage.xy_at(0, 0);
170
171 // Compute imLoc: an image locator that matches kernel locator (0,0) such that
172 // image ctrPix overlaps center of warping kernel
173 geom::Point2I subimMin = warpingOverlapBBox.getMin();
174 typename MaskedImageT::const_xy_locator const mimageLoc =
175 maskedImage.xy_at(subimMin.getX(), subimMin.getY());
176 return afw::math::convolveAtAPoint<MaskedImageT, MaskedImageT>(
177 mimageLoc, warpingKernelLoc, warpingKernelPtr->getWidth(), warpingKernelPtr->getHeight());
178}
180 afw::table::Schema &schema)
181 : _ctrl(ctrl),
182 _instFluxResultKey(
183 FluxResultKey::addFields(schema, name, "instFlux from PeakLikelihood Flux algorithm")),
184 _centroidExtractor(schema, name) {
185 _flagHandler = FlagHandler::addFields(schema, name, getFlagDefinitions());
186}
187
189 afw::image::Exposure<float> const &exposure) const {
190 // get the value from the centroid slot only
191 geom::Point2D center = _centroidExtractor(measRecord, _flagHandler);
194 MaskedImageT const &mimage = exposure.getMaskedImage();
195
206 if (!exposure.hasPsf()) {
207 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "exposure has no PSF");
208 }
209 std::shared_ptr<afw::detection::Psf const> psfPtr = exposure.getPsf();
210 if (!geom::Box2D(mimage.getBBox()).contains(center)) {
212 os << "Center = " << center << " not in exposure bbox" << mimage.getBBox();
214 }
215
216 // compute parent index and fractional offset of ctrPix: the pixel closest to "center",
217 // the centroid of the source. This assumes that the image's XY0 is (0, 0).
218 std::pair<int, double> const xCtrPixParentIndFrac = afw::image::positionToIndex(center.getX(), true);
219 std::pair<int, double> const yCtrPixParentIndFrac = afw::image::positionToIndex(center.getY(), true);
220
221 geom::Point2I ctrPixParentInd(xCtrPixParentIndFrac.first + mimage.getX0(),
222 yCtrPixParentIndFrac.first + mimage.getY0());
223 geom::Point2D ctrPixPos(afw::image::indexToPosition(ctrPixParentInd[0]),
224 afw::image::indexToPosition(ctrPixParentInd[1]));
225
226 // compute weight = 1/sum(PSF^2) for PSF at ctrPix, where PSF is normalized to a sum of 1
227 PsfAttributes psfAttr(psfPtr, ctrPixParentInd);
228 double weight = psfAttr.computeEffectiveArea();
229
230 /*
231 * Compute value of image at center of source, as shifted by a fractional pixel to center the source
232 * on ctrPix.
233 */
234 MaskedImageT::SinglePixel mimageCtrPix = computeShiftedValue(
235 mimage, _ctrl.warpingKernelName,
236 geom::Point2D(xCtrPixParentIndFrac.second, yCtrPixParentIndFrac.second), ctrPixParentInd);
237 double instFlux = mimageCtrPix.image() * weight;
238 double var = mimageCtrPix.variance() * weight * weight;
239 result.instFlux = instFlux;
240 result.instFluxErr = std::sqrt(var);
241 measRecord.set(_instFluxResultKey, result);
242 _flagHandler.setValue(measRecord, FAILURE.number, false);
243}
244
246 _flagHandler.handleFailure(measRecord, error);
247}
248
249} // namespace base
250} // namespace meas
251} // namespace lsst
int end
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
uint64_t * ptr
Definition RangeSet.cc:95
std::ostream * os
Definition Schema.cc:557
A class to contain the data, WCS, and other information needed to describe an image of the sky.
Definition Exposure.h:72
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:74
A single pixel of the same type as a MaskedImage.
Definition Pixel.h:73
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
Record class that contains measurements made on a single exposure.
Definition Source.h:78
A floating-point coordinate rectangle geometry.
Definition Box.h:413
bool contains(Point2D const &point) const noexcept
Return true if the box contains the point.
Definition Box.cc:322
An integer coordinate rectangle.
Definition Box.h:55
Point2I const getMin() const noexcept
Definition Box.h:156
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.
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.
Exception to be thrown when a measurement algorithm experiences a known failure mode.
Definition exceptions.h:48
PeakLikelihoodFluxAlgorithm(Control const &ctrl, std::string const &name, afw::table::Schema &schema)
static FlagDefinitionList const & getFlagDefinitions()
virtual void measure(afw::table::SourceRecord &measRecord, afw::image::Exposure< float > const &exposure) const
Called to measure a single child source in an image.
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.
C++ control object for peak likelihood instrument flux.
std::string warpingKernelName
"Name of warping kernel (e.g. \"lanczos4\") used to compute the peak" ;
Reports invalid arguments.
Definition Runtime.h:66
Reports when the result of an operation cannot be represented by the destination type.
Definition Runtime.h:115
T make_pair(T... args)
double indexToPosition(double ind)
Convert image index to image position.
Definition ImageUtils.h:55
int positionToIndex(double pos)
Convert image position to nearest integer index.
Definition ImageUtils.h:69
std::shared_ptr< SeparableKernel > makeWarpingKernel(std::string name)
Return a warping kernel given its name.
afw::image::MaskedImage< T >::SinglePixel computeShiftedValue(afw::image::MaskedImage< T > const &maskedImage, std::string const &warpingKernelName, geom::Point2D const &fracShift, geom::Point2I const &parentInd)
Compute the value of one pixel of an image after a fractional pixel shift Since we only want the valu...
T sqrt(T... args)
afw::table::Key< double > weight
afw::table::Key< std::string > warpingKernelName
Definition CoaddPsf.cc:354
A reusable result struct for instFlux measurements.