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
KernelCandidate.cc
Go to the documentation of this file.
1// -*- lsst-c++ -*-
11#include <stdexcept>
12
13#include "lsst/afw/math.h"
14#include "lsst/afw/image.h"
15#include "lsst/log/Log.h"
17
22
23namespace afwMath = lsst::afw::math;
24namespace afwImage = lsst::afw::image;
26
27namespace lsst {
28namespace ip {
29namespace diffim {
30
31template <typename PixelT>
32KernelCandidate<PixelT>::KernelCandidate(float const xCenter, float const yCenter,
33 MaskedImagePtr const& templateMaskedImage,
34 MaskedImagePtr const& scienceMaskedImage,
36 : lsst::afw::math::SpatialCellImageCandidate(xCenter, yCenter),
37 _templateMaskedImage(templateMaskedImage),
38 _scienceMaskedImage(scienceMaskedImage),
39 _varianceEstimate(),
40 _ps(ps.deepCopy()),
41 _source(),
42 _coreFlux(),
43 _isInitialized(false),
44 _useRegularization(false),
45 _fitForBackground(ps.getAsBool("fitForBackground")),
46 _kernelSolutionOrig(),
47 _kernelSolutionPca() {
48 /* Rank by mean core S/N in science image */
49 ImageStatistics<PixelT> imstats(ps);
50 int candidateCoreRadius = _ps->getAsInt("candidateCoreRadius");
51 try {
52 imstats.apply(*_scienceMaskedImage, candidateCoreRadius);
53 } catch (pexExcept::Exception& e) {
54 LOGL_DEBUG("TRACE2.ip.diffim.KernelCandidate",
55 "Unable to calculate core imstats for rating Candidate %d", this->getId());
56 this->setStatus(afwMath::SpatialCellCandidate::BAD);
57 return;
58 }
59
60 _coreFlux = imstats.getMean();
61 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate", "Candidate %d at %.2f %.2f with rating %.2f",
62 this->getId(), this->getXCenter(), this->getYCenter(), _coreFlux);
63}
64
65template <typename PixelT>
66KernelCandidate<PixelT>::KernelCandidate(SourcePtr const& source, MaskedImagePtr const& templateMaskedImage,
67 MaskedImagePtr const& scienceMaskedImage,
69 : lsst::afw::math::SpatialCellImageCandidate(source->getX(), source->getY()),
70 _templateMaskedImage(templateMaskedImage),
71 _scienceMaskedImage(scienceMaskedImage),
72 _varianceEstimate(),
73 _ps(ps.deepCopy()),
74 _source(source),
75 _coreFlux(source->getPsfInstFlux()),
76 _isInitialized(false),
77 _useRegularization(false),
78 _fitForBackground(ps.getAsBool("fitForBackground")),
79 _kernelSolutionOrig(),
80 _kernelSolutionPca() {
81 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate", "Candidate %d at %.2f %.2f with rating %.2f",
82 this->getId(), this->getXCenter(), this->getYCenter(), _coreFlux);
83}
84
85template <typename PixelT>
87 build(basisList, Eigen::MatrixXd());
88}
89
90template <typename PixelT>
92 Eigen::MatrixXd const& hMat) {
93 /* Examine the property set for control over the variance estimate */
95 afwImage::Image<afwImage::VariancePixel>(*(_scienceMaskedImage->getVariance()), true);
96 /* Variance estimate comes from sum of image variances */
97 var += (*(_templateMaskedImage->getVariance()));
98
99 if (_ps->getAsBool("constantVarianceWeighting")) {
100 /* Constant variance weighting */
102 float varValue;
103 if (varStats.getValue(afwMath::MEDIAN) <= 0.0)
104 varValue = 1.0;
105 else
106 varValue = varStats.getValue(afwMath::MEDIAN);
107 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate", "Candidate %d using constant variance of %.2f",
108 this->getId(), varValue);
109 var = varValue;
110 }
111
112 _varianceEstimate = VariancePtr(new afwImage::Image<afwImage::VariancePixel>(var));
113
114 try {
115 _buildKernelSolution(basisList, hMat);
116 } catch (pexExcept::Exception& e) {
117 throw e;
118 }
119
120 if (_ps->getAsBool("iterateSingleKernel") && (!(_ps->getAsBool("constantVarianceWeighting")))) {
121 afwImage::MaskedImage<PixelT> diffim = getDifferenceImage(KernelCandidate::RECENT);
122 _varianceEstimate = diffim.getVariance();
123
124 try {
125 _buildKernelSolution(basisList, hMat);
126 } catch (pexExcept::Exception& e) {
127 throw e;
128 }
129 }
130
131 _isInitialized = true;
132}
133
134template <typename PixelT>
136 Eigen::MatrixXd const& hMat) {
137 bool checkConditionNumber = _ps->getAsBool("checkConditionNumber");
138 double maxConditionNumber = _ps->getAsDouble("maxConditionNumber");
139 std::string conditionNumberType = _ps->getAsString("conditionNumberType");
141 if (conditionNumberType == "SVD") {
142 ctype = KernelSolution::SVD;
143 } else if (conditionNumberType == "EIGENVALUE") {
145 } else {
146 throw LSST_EXCEPT(pexExcept::TypeError, "conditionNumberType not recognized");
147 }
148
149 /* Do we have a regularization matrix? If so use it */
150 if (hMat.size() > 0) {
151 _useRegularization = true;
152 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate.build", "Using kernel regularization");
153
154 if (_isInitialized) {
156 new RegularizedKernelSolution<PixelT>(basisList, _fitForBackground, hMat, *_ps));
157 _kernelSolutionPca->build(*(_templateMaskedImage->getImage()), *(_scienceMaskedImage->getImage()),
158 *_varianceEstimate);
159 if (checkConditionNumber) {
160 if (_kernelSolutionPca->getConditionNumber(ctype) > maxConditionNumber) {
161 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate",
162 "Candidate %d solution has bad condition number", this->getId());
163 this->setStatus(afwMath::SpatialCellCandidate::BAD);
164 return;
165 }
166 }
167 _kernelSolutionPca->solve();
168 } else {
169 _kernelSolutionOrig = std::shared_ptr<StaticKernelSolution<PixelT> >(
170 new RegularizedKernelSolution<PixelT>(basisList, _fitForBackground, hMat, *_ps));
171 _kernelSolutionOrig->build(*(_templateMaskedImage->getImage()),
172 *(_scienceMaskedImage->getImage()), *_varianceEstimate);
173 if (checkConditionNumber) {
174 if (_kernelSolutionOrig->getConditionNumber(ctype) > maxConditionNumber) {
175 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate",
176 "Candidate %d solution has bad condition number", this->getId());
177 this->setStatus(afwMath::SpatialCellCandidate::BAD);
178 return;
179 }
180 }
181 _kernelSolutionOrig->solve();
182 }
183 } else {
184 _useRegularization = false;
185 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate.build", "Not using kernel regularization");
186 if (_isInitialized) {
188 new StaticKernelSolution<PixelT>(basisList, _fitForBackground));
189 _kernelSolutionPca->build(*(_templateMaskedImage->getImage()), *(_scienceMaskedImage->getImage()),
190 *_varianceEstimate);
191 if (checkConditionNumber) {
192 if (_kernelSolutionPca->getConditionNumber(ctype) > maxConditionNumber) {
193 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate",
194 "Candidate %d solution has bad condition number", this->getId());
195 this->setStatus(afwMath::SpatialCellCandidate::BAD);
196 return;
197 }
198 }
199 _kernelSolutionPca->solve();
200 } else {
201 _kernelSolutionOrig = std::shared_ptr<StaticKernelSolution<PixelT> >(
202 new StaticKernelSolution<PixelT>(basisList, _fitForBackground));
203 _kernelSolutionOrig->build(*(_templateMaskedImage->getImage()),
204 *(_scienceMaskedImage->getImage()), *_varianceEstimate);
205 if (checkConditionNumber) {
206 if (_kernelSolutionOrig->getConditionNumber(ctype) > maxConditionNumber) {
207 LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate",
208 "Candidate %d solution has bad condition number", this->getId());
209 this->setStatus(afwMath::SpatialCellCandidate::BAD);
210 return;
211 }
212 }
213 _kernelSolutionOrig->solve();
214 }
215 }
216}
217
218template <typename PixelT>
220 switch (cand) {
222 if (_kernelSolutionOrig)
223 return _kernelSolutionOrig->getKernel();
224 else
225 throw LSST_EXCEPT(pexExcept::RuntimeError, "Original kernel does not exist");
226 break;
228 if (_kernelSolutionPca)
229 return _kernelSolutionPca->getKernel();
230 else
231 throw LSST_EXCEPT(pexExcept::RuntimeError, "Pca kernel does not exist");
232 break;
234 if (_kernelSolutionPca)
235 return _kernelSolutionPca->getKernel();
236 else if (_kernelSolutionOrig)
237 return _kernelSolutionOrig->getKernel();
238 else
239 throw LSST_EXCEPT(pexExcept::RuntimeError, "No kernels exist");
240 break;
241 default:
242 throw std::logic_error("Invalid CandidateSwitch, cannot get kernel");
243 }
244}
245
246template <typename PixelT>
248 switch (cand) {
250 if (_kernelSolutionOrig)
251 return _kernelSolutionOrig->getBackground();
252 else
253 throw LSST_EXCEPT(pexExcept::RuntimeError, "Original kernel does not exist");
254 break;
256 if (_kernelSolutionPca)
257 return _kernelSolutionPca->getBackground();
258 else
259 throw LSST_EXCEPT(pexExcept::RuntimeError, "Pca kernel does not exist");
260 break;
262 if (_kernelSolutionPca)
263 return _kernelSolutionPca->getBackground();
264 else if (_kernelSolutionOrig)
265 return _kernelSolutionOrig->getBackground();
266 else
267 throw LSST_EXCEPT(pexExcept::RuntimeError, "No kernels exist");
268 break;
269 default:
270 throw std::logic_error("Invalid CandidateSwitch, cannot get background");
271 }
272}
273
274template <typename PixelT>
276 switch (cand) {
278 if (_kernelSolutionOrig)
279 return _kernelSolutionOrig->getKsum();
280 else
281 throw LSST_EXCEPT(pexExcept::RuntimeError, "Original kernel does not exist");
282 break;
284 if (_kernelSolutionPca)
285 return _kernelSolutionPca->getKsum();
286 else
287 throw LSST_EXCEPT(pexExcept::RuntimeError, "Pca kernel does not exist");
288 break;
290 if (_kernelSolutionPca)
291 return _kernelSolutionPca->getKsum();
292 else if (_kernelSolutionOrig)
293 return _kernelSolutionOrig->getKsum();
294 else
295 throw LSST_EXCEPT(pexExcept::RuntimeError, "No kernels exist");
296 break;
297 default:
298 throw std::logic_error("Invalid CandidateSwitch, cannot get kSum");
299 }
300}
301
302template <typename PixelT>
304 CandidateSwitch cand) const {
305 switch (cand) {
307 if (_kernelSolutionOrig)
308 return _kernelSolutionOrig->makeKernelImage();
309 else
310 throw LSST_EXCEPT(pexExcept::RuntimeError, "Original kernel does not exist");
311 break;
313 if (_kernelSolutionPca)
314 return _kernelSolutionPca->makeKernelImage();
315 else
316 throw LSST_EXCEPT(pexExcept::RuntimeError, "Pca kernel does not exist");
317 break;
319 if (_kernelSolutionPca)
320 return _kernelSolutionPca->makeKernelImage();
321 else if (_kernelSolutionOrig)
322 return _kernelSolutionOrig->makeKernelImage();
323 else
324 throw LSST_EXCEPT(pexExcept::RuntimeError, "No kernels exist");
325 break;
326 default:
327 throw std::logic_error("Invalid CandidateSwitch, cannot get kernel image");
328 }
329}
330
331template <typename PixelT>
335
336template <typename PixelT>
338 CandidateSwitch cand) const {
339 switch (cand) {
341 if (_kernelSolutionOrig)
342 return _kernelSolutionOrig;
343 else
344 throw LSST_EXCEPT(pexExcept::RuntimeError, "Original kernel does not exist");
345 break;
347 if (_kernelSolutionPca)
348 return _kernelSolutionPca;
349 else
350 throw LSST_EXCEPT(pexExcept::RuntimeError, "Pca kernel does not exist");
351 break;
353 if (_kernelSolutionPca)
354 return _kernelSolutionPca;
355 else if (_kernelSolutionOrig)
356 return _kernelSolutionOrig;
357 else
358 throw LSST_EXCEPT(pexExcept::RuntimeError, "No kernels exist");
359 break;
360 default:
361 throw std::logic_error("Invalid CandidateSwitch, cannot get solution");
362 }
363}
364
365template <typename PixelT>
367 switch (cand) {
369 if (_kernelSolutionOrig)
370 return getDifferenceImage(_kernelSolutionOrig->getKernel(),
371 _kernelSolutionOrig->getBackground());
372 else
373 throw LSST_EXCEPT(pexExcept::RuntimeError, "Original kernel does not exist");
374 break;
376 if (_kernelSolutionPca)
377 return getDifferenceImage(_kernelSolutionPca->getKernel(),
378 _kernelSolutionPca->getBackground());
379 else
380 throw LSST_EXCEPT(pexExcept::RuntimeError, "Pca kernel does not exist");
381 break;
383 if (_kernelSolutionPca)
384 return getDifferenceImage(_kernelSolutionPca->getKernel(),
385 _kernelSolutionPca->getBackground());
386 else if (_kernelSolutionOrig)
387 return getDifferenceImage(_kernelSolutionOrig->getKernel(),
388 _kernelSolutionOrig->getBackground());
389 else
390 throw LSST_EXCEPT(pexExcept::RuntimeError, "No kernels exist");
391 break;
392 default:
393 throw std::logic_error("Invalid CandidateSwitch, cannot get diffim");
394 }
395}
396
397template <typename PixelT>
399 std::shared_ptr<lsst::afw::math::Kernel> kernel, double background) {
400 /* Make diffim and set chi2 from result */
402 convolveAndSubtract(*_templateMaskedImage, *_scienceMaskedImage, *kernel, background);
403 return diffIm;
404}
405
406/***********************************************************************************************************/
407//
408// Explicit instantiations
409//
410typedef float PixelT;
411
412template class KernelCandidate<PixelT>;
413
414} // namespace diffim
415} // namespace ip
416} // namespace lsst
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
Image Subtraction helper functions.
Image Subtraction helper functions.
Class used by SpatialModelCell for spatial Kernel fitting.
Declaration of classes to store the solution for convolution kernels.
LSST DM logging module built on log4cxx.
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition Log.h:515
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
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage's variance.
float getYCenter() const
Return the object's row-centre.
Definition SpatialCell.h:91
float getXCenter() const
Return the object's column-centre.
Definition SpatialCell.h:88
int getId() const
Return the candidate's unique ID.
void setStatus(Status status)
Set the candidate's status.
A class to evaluate image statistics.
Definition Statistics.h:222
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Class for storing generic metadata.
Definition PropertySet.h:66
Class to calculate difference image statistics.
void apply(lsst::afw::image::MaskedImage< PixelT > const &image)
Class stored in SpatialCells for spatial Kernel fitting.
std::shared_ptr< ImageT const > getImage() const
afw::image::MaskedImage< PixelT > getDifferenceImage(CandidateSwitch cand)
Calculate associated difference image using internal solutions.
double getBackground(CandidateSwitch cand) const
double getKsum(CandidateSwitch cand) const
std::shared_ptr< StaticKernelSolution< PixelT > > getKernelSolution(CandidateSwitch cand) const
std::shared_ptr< afw::math::Kernel > getKernel(CandidateSwitch cand) const
Return results of kernel solution.
KernelCandidate(float const xCenter, float const yCenter, MaskedImagePtr const &templateMaskedImage, MaskedImagePtr const &scienceMaskedImage, daf::base::PropertySet const &ps)
Constructor.
std::shared_ptr< ImageT > getKernelImage(CandidateSwitch cand) const
void build(afw::math::KernelList const &basisList)
Core functionality of KernelCandidate, to build and fill a KernelSolution.
Provides consistent interface for LSST exceptions.
Definition Exception.h:107
Reports errors that are due to events beyond the control of the program.
Definition Runtime.h:104
Reports errors from accepting an object of an unexpected or inappropriate type.
Definition Runtime.h:167
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
@ MEDIAN
estimate sample median
Definition Statistics.h:60
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.