LSST Applications g070148d5b3+33e5256705,g0d53e28543+25c8b88941,g0da5cf3356+2dd1178308,g1081da9e2a+62d12e78cb,g17e5ecfddb+7e422d6136,g1c76d35bf8+ede3a706f7,g295839609d+225697d880,g2e2c1a68ba+cc1f6f037e,g2ffcdf413f+853cd4dcde,g38293774b4+62d12e78cb,g3b44f30a73+d953f1ac34,g48ccf36440+885b902d19,g4b2f1765b6+7dedbde6d2,g5320a0a9f6+0c5d6105b6,g56b687f8c9+ede3a706f7,g5c4744a4d9+ef6ac23297,g5ffd174ac0+0c5d6105b6,g6075d09f38+66af417445,g667d525e37+2ced63db88,g670421136f+2ced63db88,g71f27ac40c+2ced63db88,g774830318a+463cbe8d1f,g7876bc68e5+1d137996f1,g7985c39107+62d12e78cb,g7fdac2220c+0fd8241c05,g96f01af41f+368e6903a7,g9ca82378b8+2ced63db88,g9d27549199+ef6ac23297,gabe93b2c52+e3573e3735,gb065e2a02a+3dfbe639da,gbc3249ced9+0c5d6105b6,gbec6a3398f+0c5d6105b6,gc9534b9d65+35b9f25267,gd01420fc67+0c5d6105b6,geee7ff78d7+a14128c129,gf63283c776+ede3a706f7,gfed783d017+0c5d6105b6,w.2022.47
LSST Data Management Base Package
Loading...
Searching...
No Matches
KernelCandidate.cc
Go to the documentation of this file.
1// -*- lsst-c++ -*-
12#include "boost/timer.hpp"
13
14#include "lsst/afw/math.h"
15#include "lsst/afw/image.h"
16#include "lsst/log/Log.h"
17#include "lsst/pex/exceptions/Runtime.h"
18
19#include "lsst/ip/diffim/KernelCandidate.h"
20#include "lsst/ip/diffim/ImageSubtract.h"
21#include "lsst/ip/diffim/ImageStatistics.h"
22#include "lsst/ip/diffim/KernelSolution.h"
23
24namespace afwMath = lsst::afw::math;
25namespace afwImage = lsst::afw::image;
26namespace pexExcept = lsst::pex::exceptions;
27
28namespace lsst {
29namespace ip {
30namespace diffim {
31
32template <typename PixelT>
33KernelCandidate<PixelT>::KernelCandidate(float const xCenter, float const yCenter,
34 MaskedImagePtr const& templateMaskedImage,
35 MaskedImagePtr const& scienceMaskedImage,
36 lsst::daf::base::PropertySet const& ps)
37 : lsst::afw::math::SpatialCellImageCandidate(xCenter, yCenter),
38 _templateMaskedImage(templateMaskedImage),
39 _scienceMaskedImage(scienceMaskedImage),
40 _varianceEstimate(),
41 _ps(ps.deepCopy()),
42 _source(),
43 _coreFlux(),
44 _isInitialized(false),
45 _useRegularization(false),
46 _fitForBackground(ps.getAsBool("fitForBackground")),
47 _kernelSolutionOrig(),
48 _kernelSolutionPca() {
49 /* Rank by mean core S/N in science image */
50 ImageStatistics<PixelT> imstats(ps);
51 int candidateCoreRadius = _ps->getAsInt("candidateCoreRadius");
52 try {
53 imstats.apply(*_scienceMaskedImage, candidateCoreRadius);
54 } catch (pexExcept::Exception& e) {
55 LOGL_DEBUG("TRACE2.ip.diffim.KernelCandidate",
56 "Unable to calculate core imstats for ranking Candidate %d", this->getId());
57 this->setStatus(afwMath::SpatialCellCandidate::BAD);
58 return;
59 }
60
61 _coreFlux = imstats.getMean();
62 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate", "Candidate %d at %.2f %.2f with ranking %.2f",
63 this->getId(), this->getXCenter(), this->getYCenter(), _coreFlux);
64}
65
66template <typename PixelT>
67KernelCandidate<PixelT>::KernelCandidate(SourcePtr const& source, MaskedImagePtr const& templateMaskedImage,
68 MaskedImagePtr const& scienceMaskedImage,
69 lsst::daf::base::PropertySet const& ps)
70 : lsst::afw::math::SpatialCellImageCandidate(source->getX(), source->getY()),
71 _templateMaskedImage(templateMaskedImage),
72 _scienceMaskedImage(scienceMaskedImage),
73 _varianceEstimate(),
74 _ps(ps.deepCopy()),
75 _source(source),
76 _coreFlux(source->getPsfInstFlux()),
77 _isInitialized(false),
78 _useRegularization(false),
79 _fitForBackground(ps.getAsBool("fitForBackground")),
80 _kernelSolutionOrig(),
81 _kernelSolutionPca() {
82 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate", "Candidate %d at %.2f %.2f with ranking %.2f",
83 this->getId(), this->getXCenter(), this->getYCenter(), _coreFlux);
84}
85
86template <typename PixelT>
87void KernelCandidate<PixelT>::build(lsst::afw::math::KernelList const& basisList) {
88 build(basisList, Eigen::MatrixXd());
89}
90
91template <typename PixelT>
92void KernelCandidate<PixelT>::build(lsst::afw::math::KernelList const& basisList,
93 Eigen::MatrixXd const& hMat) {
94 /* Examine the property set for control over the variance estimate */
95 afwImage::Image<afwImage::VariancePixel> var =
96 afwImage::Image<afwImage::VariancePixel>(*(_scienceMaskedImage->getVariance()), true);
97 /* Variance estimate comes from sum of image variances */
98 var += (*(_templateMaskedImage->getVariance()));
99
100 if (_ps->getAsBool("constantVarianceWeighting")) {
101 /* Constant variance weighting */
102 afwMath::Statistics varStats = afwMath::makeStatistics(var, afwMath::MEDIAN);
103 float varValue;
104 if (varStats.getValue(afwMath::MEDIAN) <= 0.0)
105 varValue = 1.0;
106 else
107 varValue = varStats.getValue(afwMath::MEDIAN);
108 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate", "Candidate %d using constant variance of %.2f",
109 this->getId(), varValue);
110 var = varValue;
111 }
112
113 _varianceEstimate = VariancePtr(new afwImage::Image<afwImage::VariancePixel>(var));
114
115 try {
116 _buildKernelSolution(basisList, hMat);
117 } catch (pexExcept::Exception& e) {
118 throw e;
119 }
120
121 if (_ps->getAsBool("iterateSingleKernel") && (!(_ps->getAsBool("constantVarianceWeighting")))) {
122 afwImage::MaskedImage<PixelT> diffim = getDifferenceImage(KernelCandidate::RECENT);
123 _varianceEstimate = diffim.getVariance();
124
125 try {
126 _buildKernelSolution(basisList, hMat);
127 } catch (pexExcept::Exception& e) {
128 throw e;
129 }
130 }
131
132 _isInitialized = true;
133}
134
135template <typename PixelT>
136void KernelCandidate<PixelT>::_buildKernelSolution(lsst::afw::math::KernelList const& basisList,
137 Eigen::MatrixXd const& hMat) {
138 bool checkConditionNumber = _ps->getAsBool("checkConditionNumber");
139 double maxConditionNumber = _ps->getAsDouble("maxConditionNumber");
140 std::string conditionNumberType = _ps->getAsString("conditionNumberType");
141 KernelSolution::ConditionNumberType ctype;
142 if (conditionNumberType == "SVD") {
143 ctype = KernelSolution::SVD;
144 } else if (conditionNumberType == "EIGENVALUE") {
145 ctype = KernelSolution::EIGENVALUE;
146 } else {
147 throw LSST_EXCEPT(pexExcept::Exception, "conditionNumberType not recognized");
148 }
149
150 /* Do we have a regularization matrix? If so use it */
151 if (hMat.size() > 0) {
152 _useRegularization = true;
153 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate.build", "Using kernel regularization");
154
155 if (_isInitialized) {
156 _kernelSolutionPca = std::shared_ptr<StaticKernelSolution<PixelT> >(
157 new RegularizedKernelSolution<PixelT>(basisList, _fitForBackground, hMat, *_ps));
158 _kernelSolutionPca->build(*(_templateMaskedImage->getImage()), *(_scienceMaskedImage->getImage()),
159 *_varianceEstimate);
160 if (checkConditionNumber) {
161 if (_kernelSolutionPca->getConditionNumber(ctype) > maxConditionNumber) {
162 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate",
163 "Candidate %d solution has bad condition number", this->getId());
164 this->setStatus(afwMath::SpatialCellCandidate::BAD);
165 return;
166 }
167 }
168 _kernelSolutionPca->solve();
169 } else {
170 _kernelSolutionOrig = std::shared_ptr<StaticKernelSolution<PixelT> >(
171 new RegularizedKernelSolution<PixelT>(basisList, _fitForBackground, hMat, *_ps));
172 _kernelSolutionOrig->build(*(_templateMaskedImage->getImage()),
173 *(_scienceMaskedImage->getImage()), *_varianceEstimate);
174 if (checkConditionNumber) {
175 if (_kernelSolutionOrig->getConditionNumber(ctype) > maxConditionNumber) {
176 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate",
177 "Candidate %d solution has bad condition number", this->getId());
178 this->setStatus(afwMath::SpatialCellCandidate::BAD);
179 return;
180 }
181 }
182 _kernelSolutionOrig->solve();
183 }
184 } else {
185 _useRegularization = false;
186 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate.build", "Not using kernel regularization");
187 if (_isInitialized) {
188 _kernelSolutionPca = std::shared_ptr<StaticKernelSolution<PixelT> >(
189 new StaticKernelSolution<PixelT>(basisList, _fitForBackground));
190 _kernelSolutionPca->build(*(_templateMaskedImage->getImage()), *(_scienceMaskedImage->getImage()),
191 *_varianceEstimate);
192 if (checkConditionNumber) {
193 if (_kernelSolutionPca->getConditionNumber(ctype) > maxConditionNumber) {
194 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate",
195 "Candidate %d solution has bad condition number", this->getId());
196 this->setStatus(afwMath::SpatialCellCandidate::BAD);
197 return;
198 }
199 }
200 _kernelSolutionPca->solve();
201 } else {
202 _kernelSolutionOrig = std::shared_ptr<StaticKernelSolution<PixelT> >(
203 new StaticKernelSolution<PixelT>(basisList, _fitForBackground));
204 _kernelSolutionOrig->build(*(_templateMaskedImage->getImage()),
205 *(_scienceMaskedImage->getImage()), *_varianceEstimate);
206 if (checkConditionNumber) {
207 if (_kernelSolutionOrig->getConditionNumber(ctype) > maxConditionNumber) {
208 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate",
209 "Candidate %d solution has bad condition number", this->getId());
210 this->setStatus(afwMath::SpatialCellCandidate::BAD);
211 return;
212 }
213 }
214 _kernelSolutionOrig->solve();
215 }
216 }
217}
218
219template <typename PixelT>
220std::shared_ptr<lsst::afw::math::Kernel> KernelCandidate<PixelT>::getKernel(CandidateSwitch cand) const {
221 if (cand == KernelCandidate::ORIG) {
222 if (_kernelSolutionOrig)
223 return _kernelSolutionOrig->getKernel();
224 else
225 throw LSST_EXCEPT(pexExcept::Exception, "Original kernel does not exist");
226 } else if (cand == KernelCandidate::PCA) {
227 if (_kernelSolutionPca)
228 return _kernelSolutionPca->getKernel();
229 else
230 throw LSST_EXCEPT(pexExcept::Exception, "Pca kernel does not exist");
231 } else if (cand == KernelCandidate::RECENT) {
232 if (_kernelSolutionPca)
233 return _kernelSolutionPca->getKernel();
234 else if (_kernelSolutionOrig)
235 return _kernelSolutionOrig->getKernel();
236 else
237 throw LSST_EXCEPT(pexExcept::Exception, "No kernels exist");
238 } else {
239 throw LSST_EXCEPT(pexExcept::Exception, "Invalid CandidateSwitch, cannot get kernel");
240 }
241}
242
243template <typename PixelT>
244double KernelCandidate<PixelT>::getBackground(CandidateSwitch cand) const {
245 if (cand == KernelCandidate::ORIG) {
246 if (_kernelSolutionOrig)
247 return _kernelSolutionOrig->getBackground();
248 else
249 throw LSST_EXCEPT(pexExcept::Exception, "Original kernel does not exist");
250 } else if (cand == KernelCandidate::PCA) {
251 if (_kernelSolutionPca)
252 return _kernelSolutionPca->getBackground();
253 else
254 throw LSST_EXCEPT(pexExcept::Exception, "Pca kernel does not exist");
255 } else if (cand == KernelCandidate::RECENT) {
256 if (_kernelSolutionPca)
257 return _kernelSolutionPca->getBackground();
258 else if (_kernelSolutionOrig)
259 return _kernelSolutionOrig->getBackground();
260 else
261 throw LSST_EXCEPT(pexExcept::Exception, "No kernels exist");
262 } else {
263 throw LSST_EXCEPT(pexExcept::Exception, "Invalid CandidateSwitch, cannot get background");
264 }
265}
266
267template <typename PixelT>
268double KernelCandidate<PixelT>::getKsum(CandidateSwitch cand) const {
269 if (cand == KernelCandidate::ORIG) {
270 if (_kernelSolutionOrig)
271 return _kernelSolutionOrig->getKsum();
272 else
273 throw LSST_EXCEPT(pexExcept::Exception, "Original kernel does not exist");
274 } else if (cand == KernelCandidate::PCA) {
275 if (_kernelSolutionPca)
276 return _kernelSolutionPca->getKsum();
277 else
278 throw LSST_EXCEPT(pexExcept::Exception, "Pca kernel does not exist");
279 } else if (cand == KernelCandidate::RECENT) {
280 if (_kernelSolutionPca)
281 return _kernelSolutionPca->getKsum();
282 else if (_kernelSolutionOrig)
283 return _kernelSolutionOrig->getKsum();
284 else
285 throw LSST_EXCEPT(pexExcept::Exception, "No kernels exist");
286 } else {
287 throw LSST_EXCEPT(pexExcept::Exception, "Invalid CandidateSwitch, cannot get kSum");
288 }
289}
290
291template <typename PixelT>
292std::shared_ptr<typename KernelCandidate<PixelT>::ImageT> KernelCandidate<PixelT>::getKernelImage(
293 CandidateSwitch cand) const {
294 if (cand == KernelCandidate::ORIG) {
295 if (_kernelSolutionOrig)
296 return _kernelSolutionOrig->makeKernelImage();
297 else
298 throw LSST_EXCEPT(pexExcept::Exception, "Original kernel does not exist");
299 } else if (cand == KernelCandidate::PCA) {
300 if (_kernelSolutionPca)
301 return _kernelSolutionPca->makeKernelImage();
302 else
303 throw LSST_EXCEPT(pexExcept::Exception, "Pca kernel does not exist");
304 } else if (cand == KernelCandidate::RECENT) {
305 if (_kernelSolutionPca)
306 return _kernelSolutionPca->makeKernelImage();
307 else if (_kernelSolutionOrig)
308 return _kernelSolutionOrig->makeKernelImage();
309 else
310 throw LSST_EXCEPT(pexExcept::Exception, "No kernels exist");
311 } else {
312 throw LSST_EXCEPT(pexExcept::Exception, "Invalid CandidateSwitch, cannot get kernel image");
313 }
314}
315
316template <typename PixelT>
317std::shared_ptr<typename KernelCandidate<PixelT>::ImageT const> KernelCandidate<PixelT>::getImage() const {
318 return getKernelImage(KernelCandidate::ORIG);
319}
320
321template <typename PixelT>
322std::shared_ptr<StaticKernelSolution<PixelT> > KernelCandidate<PixelT>::getKernelSolution(
323 CandidateSwitch cand) const {
324 if (cand == KernelCandidate::ORIG) {
325 if (_kernelSolutionOrig)
326 return _kernelSolutionOrig;
327 else
328 throw LSST_EXCEPT(pexExcept::Exception, "Original kernel does not exist");
329 } else if (cand == KernelCandidate::PCA) {
330 if (_kernelSolutionPca)
331 return _kernelSolutionPca;
332 else
333 throw LSST_EXCEPT(pexExcept::Exception, "Pca kernel does not exist");
334 } else if (cand == KernelCandidate::RECENT) {
335 if (_kernelSolutionPca)
336 return _kernelSolutionPca;
337 else if (_kernelSolutionOrig)
338 return _kernelSolutionOrig;
339 else
340 throw LSST_EXCEPT(pexExcept::Exception, "No kernels exist");
341 } else {
342 throw LSST_EXCEPT(pexExcept::Exception, "Invalid CandidateSwitch, cannot get solution");
343 }
344}
345
346template <typename PixelT>
347lsst::afw::image::MaskedImage<PixelT> KernelCandidate<PixelT>::getDifferenceImage(CandidateSwitch cand) {
348 if (cand == KernelCandidate::ORIG) {
349 if (_kernelSolutionOrig)
350 return getDifferenceImage(_kernelSolutionOrig->getKernel(), _kernelSolutionOrig->getBackground());
351 else
352 throw LSST_EXCEPT(pexExcept::Exception, "Original kernel does not exist");
353 } else if (cand == KernelCandidate::PCA) {
354 if (_kernelSolutionPca)
355 return getDifferenceImage(_kernelSolutionPca->getKernel(), _kernelSolutionPca->getBackground());
356 else
357 throw LSST_EXCEPT(pexExcept::Exception, "Pca kernel does not exist");
358 } else if (cand == KernelCandidate::RECENT) {
359 if (_kernelSolutionPca)
360 return getDifferenceImage(_kernelSolutionPca->getKernel(), _kernelSolutionPca->getBackground());
361 else if (_kernelSolutionOrig)
362 return getDifferenceImage(_kernelSolutionOrig->getKernel(), _kernelSolutionOrig->getBackground());
363 else
364 throw LSST_EXCEPT(pexExcept::Exception, "No kernels exist");
365 } else {
366 throw LSST_EXCEPT(pexExcept::Exception, "Invalid CandidateSwitch, cannot get diffim");
367 }
368}
369
370template <typename PixelT>
371lsst::afw::image::MaskedImage<PixelT> KernelCandidate<PixelT>::getDifferenceImage(
372 std::shared_ptr<lsst::afw::math::Kernel> kernel, double background) {
373 /* Make diffim and set chi2 from result */
374 afwImage::MaskedImage<PixelT> diffIm =
375 convolveAndSubtract(*_templateMaskedImage, *_scienceMaskedImage, *kernel, background);
376 return diffIm;
377}
378
379/***********************************************************************************************************/
380//
381// Explicit instantiations
382//
383typedef float PixelT;
384
385template class KernelCandidate<PixelT>;
386
387} // namespace diffim
388} // namespace ip
389} // namespace lsst
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
ImageStatistics.h
Image Subtraction helper functions.
ImageSubtract.h
Image Subtraction helper functions.
KernelCandidate.h
Class used by SpatialModelCell for spatial Kernel fitting.
KernelSolution.h
Declaration of classes to store the solution for convolution kernels.
Log.h
LSST DM logging module built on log4cxx.
LOGL_DEBUG
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:515
lsst::afw::image::Image
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:51
lsst::afw::image::MaskedImage
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:74
lsst::afw::image::MaskedImage::getVariance
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage's variance.
Definition: MaskedImage.h:1052
lsst::afw::math::SpatialCellCandidate::getYCenter
float getYCenter() const
Return the object's row-centre.
Definition: SpatialCell.h:91
lsst::afw::math::SpatialCellCandidate::getXCenter
float getXCenter() const
Return the object's column-centre.
Definition: SpatialCell.h:88
lsst::afw::math::SpatialCellCandidate::getId
int getId() const
Return the candidate's unique ID.
Definition: SpatialCell.h:102
lsst::afw::math::SpatialCellCandidate::setStatus
void setStatus(Status status)
Set the candidate's status.
Definition: SpatialCell.cc:53
lsst::afw::math::Statistics
A class to evaluate image statistics.
Definition: Statistics.h:222
lsst::afw::math::Statistics::getValue
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1079
lsst::daf::base::PropertySet
Class for storing generic metadata.
Definition: PropertySet.h:66
lsst::ip::diffim::ImageStatistics
Class to calculate difference image statistics.
Definition: ImageStatistics.h:59
lsst::ip::diffim::ImageStatistics::apply
void apply(lsst::afw::image::MaskedImage< PixelT > const &image)
Definition: ImageStatistics.h:87
lsst::ip::diffim::KernelCandidate
Class stored in SpatialCells for spatial Kernel fitting.
Definition: KernelCandidate.h:39
lsst::ip::diffim::KernelCandidate::getImage
std::shared_ptr< ImageT const > getImage() const
Definition: KernelCandidate.cc:317
lsst::ip::diffim::KernelCandidate::getDifferenceImage
afw::image::MaskedImage< PixelT > getDifferenceImage(CandidateSwitch cand)
Calculate associated difference image using internal solutions.
Definition: KernelCandidate.cc:347
lsst::ip::diffim::KernelCandidate::getBackground
double getBackground(CandidateSwitch cand) const
Definition: KernelCandidate.cc:244
lsst::ip::diffim::KernelCandidate::getKsum
double getKsum(CandidateSwitch cand) const
Definition: KernelCandidate.cc:268
lsst::ip::diffim::KernelCandidate::getKernelSolution
std::shared_ptr< StaticKernelSolution< PixelT > > getKernelSolution(CandidateSwitch cand) const
Definition: KernelCandidate.cc:322
lsst::ip::diffim::KernelCandidate::getKernel
std::shared_ptr< afw::math::Kernel > getKernel(CandidateSwitch cand) const
Return results of kernel solution.
Definition: KernelCandidate.cc:220
lsst::ip::diffim::KernelCandidate::KernelCandidate
KernelCandidate(float const xCenter, float const yCenter, MaskedImagePtr const &templateMaskedImage, MaskedImagePtr const &scienceMaskedImage, daf::base::PropertySet const &ps)
Constructor.
Definition: KernelCandidate.cc:33
lsst::ip::diffim::KernelCandidate::getKernelImage
std::shared_ptr< ImageT > getKernelImage(CandidateSwitch cand) const
Definition: KernelCandidate.cc:292
lsst::ip::diffim::KernelCandidate::build
void build(afw::math::KernelList const &basisList)
Core functionality of KernelCandidate, to build and fill a KernelSolution.
Definition: KernelCandidate.cc:87
lsst.pex::exceptions::Exception
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
lsst::afw::math::makeStatistics
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition: Statistics.h:361
lsst::afw::math::MEDIAN
@ MEDIAN
estimate sample median
Definition: Statistics.h:60
lsst::ip::diffim::convolveAndSubtract
lsst::afw::image::MaskedImage< PixelT > convolveAndSubtract(lsst::afw::image::MaskedImage< PixelT > const &templateImage, lsst::afw::image::MaskedImage< PixelT > const &scienceMaskedImage, lsst::afw::math::Kernel const &convolutionKernel, BackgroundT background, bool invert=true)
Execute fundamental task of convolving template and subtracting it from science image.
Definition: ImageSubtract.cc:115
Generated on Thu Nov 24 2022 08:01:22 for LSST Applications by doxygen 1.9.5