LSSTApplications  20.0.0
LSSTDataManagementBasePackage
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 
26 #include "lsst/afw/detection/Psf.h"
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 
34 #include "lsst/afw/table/Source.h"
35 
36 namespace lsst {
37 namespace meas {
38 namespace base {
39 namespace {
40 FlagDefinitionList flagDefinitions;
41 } // namespace
42 
43 FlagDefinition const PeakLikelihoodFluxAlgorithm::FAILURE = flagDefinitions.addFailureFlag();
44 
46 
47 namespace {
48 
49 /************************************************************************************************************/
66 class PsfAttributes {
67 public:
68  enum Method {
69  ADAPTIVE_MOMENT,
70  FIRST_MOMENT,
71  SECOND_MOMENT,
72  NOISE_EQUIVALENT,
73  BICKERTON
74  };
75 
76  PsfAttributes(CONST_PTR(afw::detection::Psf) psf, int const iX, int const iY);
77  PsfAttributes(CONST_PTR(afw::detection::Psf) psf, geom::Point2I const &cen);
78 
79  double computeGaussianWidth(Method how = ADAPTIVE_MOMENT) const;
80  double computeEffectiveArea() const;
81 
82 private:
83  PTR(afw::image::Image<double>) _psfImage;
84 };
85 
89 PsfAttributes::PsfAttributes(CONST_PTR(afw::detection::Psf) 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 
100 PsfAttributes::PsfAttributes(
101  CONST_PTR(afw::detection::Psf) psf,
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 
111 double 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 
133 template <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 }
181  : _ctrl(ctrl),
182  _instFluxResultKey(
183  FluxResultKey::addFields(schema, name, "instFlux from PeakLikelihood Flux algorithm")),
184  _centroidExtractor(schema, name) {
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);
193  typedef afw::image::Exposure<float>::MaskedImageT MaskedImageT;
194  MaskedImageT const &mimage = exposure.getMaskedImage();
195 
206  if (!exposure.hasPsf()) {
207  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "exposure has no PSF");
208  }
209  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
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, yCtrPixParentIndFrac.first);
222  geom::Point2D ctrPixPos(afw::image::indexToPosition(ctrPixParentInd[0]),
223  afw::image::indexToPosition(ctrPixParentInd[1]));
224 
225  // compute weight = 1/sum(PSF^2) for PSF at ctrPix, where PSF is normalized to a sum of 1
226  PsfAttributes psfAttr(psfPtr, ctrPixParentInd);
227  double weight = psfAttr.computeEffectiveArea();
228 
229  /*
230  * Compute value of image at center of source, as shifted by a fractional pixel to center the source
231  * on ctrPix.
232  */
233  MaskedImageT::SinglePixel mimageCtrPix = computeShiftedValue(
234  mimage, _ctrl.warpingKernelName,
235  geom::Point2D(xCtrPixParentIndFrac.second, yCtrPixParentIndFrac.second), ctrPixParentInd);
236  double instFlux = mimageCtrPix.image() * weight;
237  double var = mimageCtrPix.variance() * weight * weight;
238  result.instFlux = instFlux;
239  result.instFluxErr = std::sqrt(var);
240  measRecord.set(_instFluxResultKey, result);
241  _flagHandler.setValue(measRecord, FAILURE.number, false);
242 }
243 
245  _flagHandler.handleFailure(measRecord, error);
246 }
247 
248 } // namespace base
249 } // namespace meas
250 } // namespace lsst
schema
table::Schema schema
Definition: Amplifier.cc:115
lsst::meas::base::PeakLikelihoodFluxAlgorithm::FAILURE
static FlagDefinition const FAILURE
Definition: PeakLikelihoodFlux.h:75
lsst::afw::image::Exposure::getPsf
std::shared_ptr< lsst::afw::detection::Psf const > getPsf() const
Return the Exposure's Psf object.
Definition: Exposure.h:307
CONST_PTR
#define CONST_PTR(...)
A shared pointer to a const object.
Definition: base.h:47
Box.h
lsst::meas::base::FlagDefinition::number
std::size_t number
Definition: FlagHandler.h:54
lsst::afw::image::MaskedImage::xy_at
xy_locator xy_at(int x, int y) const
Return an xy_locator at the point (x, y)
Definition: MaskedImage.h:1253
std::string
STL class.
Psf.h
lsst::meas::base::FluxResult
A reusable result struct for instFlux measurements.
Definition: FluxUtilities.h:41
lsst::log.log.logContinued.error
def error(fmt, *args)
Definition: logContinued.py:210
lsst::afw::image::positionToIndex
int positionToIndex(double pos)
Convert image position to nearest integer index.
Definition: ImageUtils.h:69
lsst::meas::base::FlagHandler::handleFailure
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
lsst::afw::table::SourceRecord
Record class that contains measurements made on a single exposure.
Definition: Source.h:80
lsst::afw::image::pixel::SinglePixel
A single pixel of the same type as a MaskedImage.
Definition: Pixel.h:73
lsst::afw::image::Exposure< float >
lsst::meas::base::PeakLikelihoodFluxControl
C++ control object for peak likelihood instrument flux.
Definition: PeakLikelihoodFlux.h:59
base
Definition: __init__.py:1
std::pair
lsst::geom::Box2D::contains
bool contains(Point2D const &point) const noexcept
Return true if the box contains the point.
Definition: Box.cc:322
lsst::pex::exceptions::RangeError
Reports when the result of an operation cannot be represented by the destination type.
Definition: Runtime.h:115
psf
Key< int > psf
Definition: Exposure.cc:65
lsst::sphgeom::abs
Angle abs(Angle const &a)
Definition: Angle.h:106
lsst::meas::base::PeakLikelihoodFluxAlgorithm::fail
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: PeakLikelihoodFlux.cc:244
lsst::afw::math::makeWarpingKernel
std::shared_ptr< SeparableKernel > makeWarpingKernel(std::string name)
Return a warping kernel given its name.
Definition: warpExposure.cc:161
lsst::meas::base::MeasurementError
Exception to be thrown when a measurement algorithm experiences a known failure mode.
Definition: exceptions.h:48
lsst::afw::image::Exposure::hasPsf
bool hasPsf() const
Does this Exposure have a Psf?
Definition: Exposure.h:310
lsst::afw::table::Schema
Defines the fields and offsets for a table.
Definition: Schema.h:50
lsst::afw::image::MaskedImage::getBBox
lsst::geom::Box2I getBBox(ImageOrigin const origin=PARENT) const
Definition: MaskedImage.h:1097
lsst::geom::Box2I::getMin
Point2I const getMin() const noexcept
Definition: Box.h:156
lsst::afw::image::indexToPosition
double indexToPosition(double ind)
Convert image index to image position.
Definition: ImageUtils.h:55
geom.h
lsst::afw::geom.transform.transformContinued.name
string name
Definition: transformContinued.py:32
lsst::meas::base::PeakLikelihoodFluxAlgorithm::PeakLikelihoodFluxAlgorithm
PeakLikelihoodFluxAlgorithm(Control const &ctrl, std::string const &name, afw::table::Schema &schema)
Definition: PeakLikelihoodFlux.cc:179
lsst::meas::base::FluxResultKey
A FunctorKey for FluxResult.
Definition: FluxUtilities.h:59
lsst::afw::image::Exposure::getMaskedImage
MaskedImageT getMaskedImage()
Return the MaskedImage.
Definition: Exposure.h:228
std::sqrt
T sqrt(T... args)
warpingKernelName
afw::table::Key< std::string > warpingKernelName
Definition: CoaddPsf.cc:341
end
int end
Definition: BoundedField.cc:105
lsst::meas::base::FlagDefinitionList
vector-type utility class to build a collection of FlagDefinitions
Definition: FlagHandler.h:60
math.h
lsst::meas::base::PeakLikelihoodFluxControl::warpingKernelName
std::string warpingKernelName
"Name of warping kernel (e.g. \"lanczos4") used to compute the peak" ;
Definition: PeakLikelihoodFlux.h:62
lsst::meas::base::FlagHandler::setValue
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
lsst::afw::image::MaskedImage
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
lsst::geom::Box2I::contains
bool contains(Point2I const &point) const noexcept
Return true if the box contains the point.
Definition: Box.cc:114
lsst::meas::base::PeakLikelihoodFluxAlgorithm::measure
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: PeakLikelihoodFlux.cc:188
lsst::meas::base::computeShiftedValue
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...
Definition: PeakLikelihoodFlux.cc:134
ptr
uint64_t * ptr
Definition: RangeSet.cc:88
image.h
result
py::object result
Definition: _schema.cc:429
weight
afw::table::Key< double > weight
Definition: CoaddBoundedField.cc:166
lsst::afw::image::ImageBase::x_iterator
_view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: ImageBase.h:133
Source.h
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
std::ostringstream
STL class.
lsst::geom
Definition: geomOperators.dox:4
os
std::ostream * os
Definition: Schema.cc:746
lsst::pex::exceptions::InvalidParameterError
Reports invalid arguments.
Definition: Runtime.h:66
PTR
#define PTR(...)
Definition: base.h:41
lsst::geom::Point< int, 2 >
ast::PointD
std::vector< double > PointD
Vector of double; used for bounds, points.
Definition: base.h:57
lsst::geom::Box2I
An integer coordinate rectangle.
Definition: Box.h:55
std::make_pair
T make_pair(T... args)
lsst::afw::detection::Psf
A polymorphic base class for representing an image's Point Spread Function.
Definition: Psf.h:76
lsst::afw::table::BaseRecord::set
void set(Key< T > const &key, U const &value)
Set value of a field for the given key.
Definition: BaseRecord.h:164
PeakLikelihoodFlux.h
lsst::geom::Box2D
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
lsst::afw::image::Image
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:58
lsst::geom::Extent< int, 2 >
lsst::afw::math::SeparableKernel
A kernel described by a pair of functions: func(x, y) = colFunc(x) * rowFunc(y)
Definition: Kernel.h:907
lsst::meas::base::FlagHandler::addFields
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
exceptions.h
lsst::meas::base::PeakLikelihoodFluxAlgorithm::getFlagDefinitions
static FlagDefinitionList const & getFlagDefinitions()
Definition: PeakLikelihoodFlux.cc:45