LSSTApplications  20.0.0
LSSTDataManagementBasePackage
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 
24 namespace afwMath = lsst::afw::math;
25 namespace afwImage = lsst::afw::image;
27 
28 namespace lsst {
29 namespace ip {
30 namespace diffim {
31 
32 template <typename PixelT>
33 KernelCandidate<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 
66 template <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 
86 template <typename PixelT>
88  build(basisList, Eigen::MatrixXd());
89 }
90 
91 template <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 
135 template <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) {
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 
219 template <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 
243 template <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 
267 template <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 
291 template <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 
316 template <typename PixelT>
318  return getKernelImage(KernelCandidate::ORIG);
319 }
320 
321 template <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 
346 template <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 
370 template <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 //
383 typedef float PixelT;
384 
385 template class KernelCandidate<PixelT>;
386 
387 } // namespace diffim
388 } // namespace ip
389 } // namespace lsst
lsst::afw::image
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Definition: imageAlgorithm.dox:1
lsst::ip::diffim::ImageStatistics::getMean
double getMean() const
Definition: ImageStatistics.h:130
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
std::string
STL class.
std::shared_ptr< afw::image::MaskedImage< PixelT > >
lsst::ip::diffim::KernelSolution::SVD
@ SVD
Definition: KernelSolution.h:47
lsst::afw::math::SpatialCellCandidate::setStatus
void setStatus(Status status)
Set the candidate's status.
Definition: SpatialCell.cc:54
lsst::ip::diffim::KernelCandidate::PCA
@ PCA
Definition: KernelCandidate.h:51
lsst::ip::diffim::ImageStatistics
Class to calculate difference image statistics.
Definition: ImageStatistics.h:59
std::vector< std::shared_ptr< Kernel > >
lsst::ip::diffim::KernelCandidate
Class stored in SpatialCells for spatial Kernel fitting.
Definition: KernelCandidate.h:39
lsst::ip::diffim::PixelT
float PixelT
Definition: KernelCandidate.cc:383
lsst::afw
Definition: imageAlgorithm.dox:1
KernelSolution.h
Declaration of classes to store the solution for convolution kernels.
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
lsst::ip::diffim::KernelCandidate::CandidateSwitch
CandidateSwitch
Definition: KernelCandidate.h:49
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:1056
lsst::ip::diffim::KernelCandidate::getDifferenceImage
afw::image::MaskedImage< PixelT > getDifferenceImage(CandidateSwitch cand)
Calculate associated difference image using internal solutions.
Definition: KernelCandidate.cc:347
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:354
math.h
lsst::afw::math::SpatialCellCandidate::BAD
@ BAD
Definition: SpatialCell.h:72
lsst::afw::math::SpatialCellCandidate::getXCenter
float getXCenter() const
Return the object's column-centre.
Definition: SpatialCell.h:90
lsst::afw::image::MaskedImage< PixelT >
ast::detail::source
const char * source()
Source function that allows astChannel to source from a Stream.
Definition: Stream.h:224
lsst::ip::diffim::KernelSolution::ConditionNumberType
ConditionNumberType
Definition: KernelSolution.h:45
lsst.pipe.drivers.visualizeVisit.background
background
Definition: visualizeVisit.py:37
image
afw::table::Key< afw::table::Array< ImagePixelT > > image
Definition: HeavyFootprint.cc:216
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::ORIG
@ ORIG
Definition: KernelCandidate.h:50
lsst::afw::math::SpatialCellCandidate::getId
int getId() const
Return the candidate's unique ID.
Definition: SpatialCell.h:104
image.h
ImageSubtract.h
Image Subtraction helper functions.
lsst::ip::diffim::KernelCandidate::getImage
boost::shared_ptr< ImageT const > getImage() const
Definition: KernelCandidate.cc:317
lsst::ip::diffim::KernelCandidate::getKernelImage
boost::shared_ptr< ImageT > getKernelImage(CandidateSwitch cand) const
Definition: KernelCandidate.cc:292
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
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::afw::image::MaskedImage::getVariance
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage's variance.
Definition: MaskedImage.h:1090
lsst::ip::diffim::ImageStatistics::apply
void apply(lsst::afw::image::MaskedImage< PixelT > const &image)
Definition: ImageStatistics.h:87
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
lsst::ip::diffim::KernelCandidate::getKernelSolution
std::shared_ptr< StaticKernelSolution< PixelT > > getKernelSolution(CandidateSwitch cand) const
Definition: KernelCandidate.cc:322
Runtime.h
lsst::ip::diffim::KernelCandidate::getBackground
double getBackground(CandidateSwitch cand) const
Definition: KernelCandidate.cc:244
LOGL_DEBUG
#define LOGL_DEBUG(logger, message...)
Definition: Log.h:504
lsst::afw::math
Definition: statistics.dox:6
KernelCandidate.h
Class used by SpatialModelCell for spatial Kernel fitting.
lsst::daf::base::PropertySet
Class for storing generic metadata.
Definition: PropertySet.h:67
lsst::ip::diffim::KernelCandidate::RECENT
@ RECENT
Definition: KernelCandidate.h:52
lsst::afw::math::Statistics
Definition: Statistics.h:215
lsst::pex::exceptions
Definition: Exception.h:37
lsst::pex::exceptions::Exception
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
lsst::ip::diffim::KernelSolution::EIGENVALUE
@ EIGENVALUE
Definition: KernelSolution.h:46
lsst::afw::image::Image
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:58
lsst::afw::math::SpatialCellCandidate::getYCenter
float getYCenter() const
Return the object's row-centre.
Definition: SpatialCell.h:93
Log.h
LSST DM logging module built on log4cxx.
lsst::afw::math::MEDIAN
@ MEDIAN
estimate sample median
Definition: Statistics.h:70
lsst::ip::diffim::KernelCandidate::getKsum
double getKsum(CandidateSwitch cand) const
Definition: KernelCandidate.cc:268
ImageStatistics.h
Image Subtraction helper functions.