LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
LSST Data Management Base Package
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"
18
23
24namespace afwMath = lsst::afw::math;
25namespace afwImage = lsst::afw::image;
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,
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());
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>
68 MaskedImagePtr const& scienceMaskedImage,
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>
88 build(basisList, Eigen::MatrixXd());
89}
90
91template <typename PixelT>
93 Eigen::MatrixXd const& hMat) {
94 /* Examine the property set for control over the variance estimate */
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 */
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>
137 Eigen::MatrixXd const& hMat) {
138 bool checkConditionNumber = _ps->getAsBool("checkConditionNumber");
139 double maxConditionNumber = _ps->getAsDouble("maxConditionNumber");
140 std::string conditionNumberType = _ps->getAsString("conditionNumberType");
142 if (conditionNumberType == "SVD") {
143 ctype = KernelSolution::SVD;
144 } else if (conditionNumberType == "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) {
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) {
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>
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>
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>
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>
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>
318 return getKernelImage(KernelCandidate::ORIG);
319}
320
321template <typename PixelT>
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>
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>
373 /* Make diffim and set chi2 from result */
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
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
afw::table::Key< afw::table::Array< ImagePixelT > > image
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:73
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage's variance.
Definition: MaskedImage.h:1051
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.
Definition: SpatialCell.h:102
void setStatus(Status status)
Set the candidate's status.
Definition: SpatialCell.cc:53
A class to evaluate image statistics.
Definition: Statistics.h:220
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1047
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
const char * source()
Source function that allows astChannel to source from a Stream.
Definition: Stream.h:224
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
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:359
@ MEDIAN
estimate sample median
Definition: Statistics.h:69
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.
A base class for image defects.