LSSTApplications  11.0-13-gbb96280,12.1+18,12.1+7,12.1-1-g14f38d3+72,12.1-1-g16c0db7+5,12.1-1-g5961e7a+84,12.1-1-ge22e12b+23,12.1-11-g06625e2+4,12.1-11-g0d7f63b+4,12.1-19-gd507bfc,12.1-2-g7dda0ab+38,12.1-2-gc0bc6ab+81,12.1-21-g6ffe579+2,12.1-21-gbdb6c2a+4,12.1-24-g941c398+5,12.1-3-g57f6835+7,12.1-3-gf0736f3,12.1-37-g3ddd237,12.1-4-gf46015e+5,12.1-5-g06c326c+20,12.1-5-g648ee80+3,12.1-5-gc2189d7+4,12.1-6-ga608fc0+1,12.1-7-g3349e2a+5,12.1-7-gfd75620+9,12.1-9-g577b946+5,12.1-9-gc4df26a+10
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;
26 namespace pexExcept = lsst::pex::exceptions;
27 
28 namespace lsst {
29 namespace ip {
30 namespace diffim {
31 
32  template <typename PixelT>
34  float const xCenter,
35  float const yCenter,
36  MaskedImagePtr const& templateMaskedImage,
37  MaskedImagePtr const& scienceMaskedImage,
38  lsst::pex::policy::Policy const& policy
39  ) :
40  lsst::afw::math::SpatialCellImageCandidate(xCenter, yCenter),
41  _templateMaskedImage(templateMaskedImage),
42  _scienceMaskedImage(scienceMaskedImage),
43  _varianceEstimate(),
44  _policy(policy),
45  _source(),
46  _coreFlux(),
47  _isInitialized(false),
48  _useRegularization(false),
49  _fitForBackground(_policy.getBool("fitForBackground")),
50  _kernelSolutionOrig(),
51  _kernelSolutionPca()
52  {
53 
54  /* Rank by mean core S/N in science image */
56  int candidateCoreRadius = _policy.getInt("candidateCoreRadius");
57  try {
58  imstats.apply(*_scienceMaskedImage, candidateCoreRadius);
59  } catch (pexExcept::Exception& e) {
60  LOGL_DEBUG("TRACE2.ip.diffim.KernelCandidate",
61  "Unable to calculate core imstats for ranking Candidate %d", this->getId());
63  return;
64  }
65 
66  _coreFlux = imstats.getMean();
67  LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate",
68  "Candidate %d at %.2f %.2f with ranking %.2f",
69  this->getId(), this->getXCenter(), this->getYCenter(), _coreFlux);
70  }
71 
72  template <typename PixelT>
74  SourcePtr const& source,
75  MaskedImagePtr const& templateMaskedImage,
76  MaskedImagePtr const& scienceMaskedImage,
77  lsst::pex::policy::Policy const& policy
78  ) :
79  lsst::afw::math::SpatialCellImageCandidate(source->getX(), source->getY()),
80  _templateMaskedImage(templateMaskedImage),
81  _scienceMaskedImage(scienceMaskedImage),
82  _varianceEstimate(),
83  _policy(policy),
84  _source(source),
85  _coreFlux(source->getPsfFlux()),
86  _isInitialized(false),
87  _useRegularization(false),
88  _fitForBackground(_policy.getBool("fitForBackground")),
89  _kernelSolutionOrig(),
90  _kernelSolutionPca()
91  {
92  LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate",
93  "Candidate %d at %.2f %.2f with ranking %.2f",
94  this->getId(), this->getXCenter(), this->getYCenter(), _coreFlux);
95  }
96 
97  template <typename PixelT>
99  lsst::afw::math::KernelList const& basisList
100  ) {
101  build(basisList, std::shared_ptr<Eigen::MatrixXd>());
102  }
103 
104  template <typename PixelT>
106  lsst::afw::math::KernelList const& basisList,
107  std::shared_ptr<Eigen::MatrixXd> hMat
108  ) {
109 
110  /* Examine the policy for control over the variance estimate */
112  afwImage::Image<afwImage::VariancePixel>(*(_scienceMaskedImage->getVariance()), true);
113  /* Variance estimate comes from sum of image variances */
114  var += (*(_templateMaskedImage->getVariance()));
115 
116  if (_policy.getBool("constantVarianceWeighting")) {
117  /* Constant variance weighting */
119  float varValue;
120  if (varStats.getValue(afwMath::MEDIAN) <= 0.0)
121  varValue = 1.0;
122  else
123  varValue = varStats.getValue(afwMath::MEDIAN);
124  LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate",
125  "Candidate %d using constant variance of %.2f", this->getId(), varValue);
126  var = varValue;
127 
128  }
129 
130  _varianceEstimate = VariancePtr( new afwImage::Image<afwImage::VariancePixel>(var) );
131 
132  try {
133  _buildKernelSolution(basisList, hMat);
134  } catch (pexExcept::Exception &e) {
135  throw e;
136  }
137 
138  if (_policy.getBool("iterateSingleKernel") && (!(_policy.getBool("constantVarianceWeighting")))) {
139  afwImage::MaskedImage<PixelT> diffim = getDifferenceImage(KernelCandidate::RECENT);
140  _varianceEstimate = diffim.getVariance();
141 
142  try {
143  _buildKernelSolution(basisList, hMat);
144  } catch (pexExcept::Exception &e) {
145  throw e;
146  }
147  }
148 
149  _isInitialized = true;
150 
151  }
152 
153  template <typename PixelT>
155  std::shared_ptr<Eigen::MatrixXd> hMat)
156  {
157  bool checkConditionNumber = _policy.getBool("checkConditionNumber");
158  double maxConditionNumber = _policy.getDouble("maxConditionNumber");
159  std::string conditionNumberType = _policy.getString("conditionNumberType");
161  if (conditionNumberType == "SVD") {
162  ctype = KernelSolution::SVD;
163  }
164  else if (conditionNumberType == "EIGENVALUE") {
166  }
167  else {
168  throw LSST_EXCEPT(pexExcept::Exception, "conditionNumberType not recognized");
169  }
170 
171  /* Do we have a regularization matrix? If so use it */
172  if (hMat) {
173  _useRegularization = true;
174  LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate.build",
175  "Using kernel regularization");
176 
177  if (_isInitialized) {
178  _kernelSolutionPca = std::shared_ptr<StaticKernelSolution<PixelT> >(
179  new RegularizedKernelSolution<PixelT>(basisList, _fitForBackground, hMat, _policy)
180  );
181  _kernelSolutionPca->build(*(_templateMaskedImage->getImage()),
182  *(_scienceMaskedImage->getImage()),
183  *_varianceEstimate);
184  if (checkConditionNumber) {
185  if (_kernelSolutionPca->getConditionNumber(ctype) > maxConditionNumber) {
186  LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate",
187  "Candidate %d solution has bad condition number",
188  this->getId());
189  this->setStatus(afwMath::SpatialCellCandidate::BAD);
190  return;
191  }
192  }
193  _kernelSolutionPca->solve();
194  }
195  else {
196  _kernelSolutionOrig = std::shared_ptr<StaticKernelSolution<PixelT> >(
197  new RegularizedKernelSolution<PixelT>(basisList, _fitForBackground, hMat, _policy)
198  );
199  _kernelSolutionOrig->build(*(_templateMaskedImage->getImage()),
200  *(_scienceMaskedImage->getImage()),
201  *_varianceEstimate);
202  if (checkConditionNumber) {
203  if (_kernelSolutionOrig->getConditionNumber(ctype) > maxConditionNumber) {
204  LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate",
205  "Candidate %d solution has bad condition number",
206  this->getId());
207  this->setStatus(afwMath::SpatialCellCandidate::BAD);
208  return;
209  }
210  }
211  _kernelSolutionOrig->solve();
212  }
213  }
214  else {
215  _useRegularization = false;
216  LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate.build",
217  "Not using kernel regularization");
218  if (_isInitialized) {
219  _kernelSolutionPca = std::shared_ptr<StaticKernelSolution<PixelT> >(
220  new StaticKernelSolution<PixelT>(basisList, _fitForBackground)
221  );
222  _kernelSolutionPca->build(*(_templateMaskedImage->getImage()),
223  *(_scienceMaskedImage->getImage()),
224  *_varianceEstimate);
225  if (checkConditionNumber) {
226  if (_kernelSolutionPca->getConditionNumber(ctype) > maxConditionNumber) {
227  LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate",
228  "Candidate %d solution has bad condition number",
229  this->getId());
230  this->setStatus(afwMath::SpatialCellCandidate::BAD);
231  return;
232  }
233  }
234  _kernelSolutionPca->solve();
235  }
236  else {
237  _kernelSolutionOrig = std::shared_ptr<StaticKernelSolution<PixelT> >(
238  new StaticKernelSolution<PixelT>(basisList, _fitForBackground)
239  );
240  _kernelSolutionOrig->build(*(_templateMaskedImage->getImage()),
241  *(_scienceMaskedImage->getImage()),
242  *_varianceEstimate);
243  if (checkConditionNumber) {
244  if (_kernelSolutionOrig->getConditionNumber(ctype) > maxConditionNumber) {
245  LOGL_DEBUG("TRACE4.ip.diffim.KernelCandidate",
246  "Candidate %d solution has bad condition number",
247  this->getId());
248  this->setStatus(afwMath::SpatialCellCandidate::BAD);
249  return;
250  }
251  }
252  _kernelSolutionOrig->solve();
253  }
254  }
255  }
256 
257  template <typename PixelT>
259  if (cand == KernelCandidate::ORIG) {
260  if (_kernelSolutionOrig)
261  return _kernelSolutionOrig->getKernel();
262  else
263  throw LSST_EXCEPT(pexExcept::Exception, "Original kernel does not exist");
264  }
265  else if (cand == KernelCandidate::PCA) {
266  if (_kernelSolutionPca)
267  return _kernelSolutionPca->getKernel();
268  else
269  throw LSST_EXCEPT(pexExcept::Exception, "Pca kernel does not exist");
270  }
271  else if (cand == KernelCandidate::RECENT) {
272  if (_kernelSolutionPca)
273  return _kernelSolutionPca->getKernel();
274  else if (_kernelSolutionOrig)
275  return _kernelSolutionOrig->getKernel();
276  else
277  throw LSST_EXCEPT(pexExcept::Exception, "No kernels exist");
278  }
279  else {
280  throw LSST_EXCEPT(pexExcept::Exception, "Invalid CandidateSwitch, cannot get kernel");
281  }
282  }
283 
284  template <typename PixelT>
286  if (cand == KernelCandidate::ORIG) {
287  if (_kernelSolutionOrig)
288  return _kernelSolutionOrig->getBackground();
289  else
290  throw LSST_EXCEPT(pexExcept::Exception, "Original kernel does not exist");
291  }
292  else if (cand == KernelCandidate::PCA) {
293  if (_kernelSolutionPca)
294  return _kernelSolutionPca->getBackground();
295  else
296  throw LSST_EXCEPT(pexExcept::Exception, "Pca kernel does not exist");
297  }
298  else if (cand == KernelCandidate::RECENT) {
299  if (_kernelSolutionPca)
300  return _kernelSolutionPca->getBackground();
301  else if (_kernelSolutionOrig)
302  return _kernelSolutionOrig->getBackground();
303  else
304  throw LSST_EXCEPT(pexExcept::Exception, "No kernels exist");
305  }
306  else {
307  throw LSST_EXCEPT(pexExcept::Exception, "Invalid CandidateSwitch, cannot get background");
308  }
309  }
310 
311  template <typename PixelT>
313  if (cand == KernelCandidate::ORIG) {
314  if (_kernelSolutionOrig)
315  return _kernelSolutionOrig->getKsum();
316  else
317  throw LSST_EXCEPT(pexExcept::Exception, "Original kernel does not exist");
318  }
319  else if (cand == KernelCandidate::PCA) {
320  if (_kernelSolutionPca)
321  return _kernelSolutionPca->getKsum();
322  else
323  throw LSST_EXCEPT(pexExcept::Exception, "Pca kernel does not exist");
324  }
325  else if (cand == KernelCandidate::RECENT) {
326  if (_kernelSolutionPca)
327  return _kernelSolutionPca->getKsum();
328  else if (_kernelSolutionOrig)
329  return _kernelSolutionOrig->getKsum();
330  else
331  throw LSST_EXCEPT(pexExcept::Exception, "No kernels exist");
332  }
333  else {
334  throw LSST_EXCEPT(pexExcept::Exception, "Invalid CandidateSwitch, cannot get kSum");
335  }
336  }
337 
338  template <typename PixelT>
340  CandidateSwitch cand) const {
341  if (cand == KernelCandidate::ORIG) {
342  if (_kernelSolutionOrig)
343  return _kernelSolutionOrig->makeKernelImage();
344  else
345  throw LSST_EXCEPT(pexExcept::Exception, "Original kernel does not exist");
346  }
347  else if (cand == KernelCandidate::PCA) {
348  if (_kernelSolutionPca)
349  return _kernelSolutionPca->makeKernelImage();
350  else
351  throw LSST_EXCEPT(pexExcept::Exception, "Pca kernel does not exist");
352  }
353  else if (cand == KernelCandidate::RECENT) {
354  if (_kernelSolutionPca)
355  return _kernelSolutionPca->makeKernelImage();
356  else if (_kernelSolutionOrig)
357  return _kernelSolutionOrig->makeKernelImage();
358  else
359  throw LSST_EXCEPT(pexExcept::Exception, "No kernels exist");
360  }
361  else {
362  throw LSST_EXCEPT(pexExcept::Exception, "Invalid CandidateSwitch, cannot get kernel image");
363  }
364  }
365 
366  template <typename PixelT>
368  return getKernelImage(KernelCandidate::ORIG);
369  }
370 
371  template <typename PixelT>
372  std::shared_ptr<StaticKernelSolution<PixelT> > KernelCandidate<PixelT>::getKernelSolution(
373  CandidateSwitch cand) const {
374  if (cand == KernelCandidate::ORIG) {
375  if (_kernelSolutionOrig)
376  return _kernelSolutionOrig;
377  else
378  throw LSST_EXCEPT(pexExcept::Exception, "Original kernel does not exist");
379  }
380  else if (cand == KernelCandidate::PCA) {
381  if (_kernelSolutionPca)
382  return _kernelSolutionPca;
383  else
384  throw LSST_EXCEPT(pexExcept::Exception, "Pca kernel does not exist");
385  }
386  else if (cand == KernelCandidate::RECENT) {
387  if (_kernelSolutionPca)
388  return _kernelSolutionPca;
389  else if (_kernelSolutionOrig)
390  return _kernelSolutionOrig;
391  else
392  throw LSST_EXCEPT(pexExcept::Exception, "No kernels exist");
393  }
394  else {
395  throw LSST_EXCEPT(pexExcept::Exception, "Invalid CandidateSwitch, cannot get solution");
396  }
397  }
398 
399  template <typename PixelT>
401  CandidateSwitch cand) {
402  if (cand == KernelCandidate::ORIG) {
403  if (_kernelSolutionOrig)
404  return getDifferenceImage(_kernelSolutionOrig->getKernel(),
405  _kernelSolutionOrig->getBackground());
406  else
407  throw LSST_EXCEPT(pexExcept::Exception, "Original kernel does not exist");
408  }
409  else if (cand == KernelCandidate::PCA) {
410  if (_kernelSolutionPca)
411  return getDifferenceImage(_kernelSolutionPca->getKernel(),
412  _kernelSolutionPca->getBackground());
413  else
414  throw LSST_EXCEPT(pexExcept::Exception, "Pca kernel does not exist");
415  }
416  else if (cand == KernelCandidate::RECENT) {
417  if (_kernelSolutionPca)
418  return getDifferenceImage(_kernelSolutionPca->getKernel(),
419  _kernelSolutionPca->getBackground());
420  else if (_kernelSolutionOrig)
421  return getDifferenceImage(_kernelSolutionOrig->getKernel(),
422  _kernelSolutionOrig->getBackground());
423  else
424  throw LSST_EXCEPT(pexExcept::Exception, "No kernels exist");
425  }
426  else {
427  throw LSST_EXCEPT(pexExcept::Exception, "Invalid CandidateSwitch, cannot get diffim");
428  }
429  }
430 
431  template <typename PixelT>
434  double background
435  ) {
436  /* Make diffim and set chi2 from result */
437  afwImage::MaskedImage<PixelT> diffIm = convolveAndSubtract(*_templateMaskedImage,
438  *_scienceMaskedImage,
439  *kernel,
440  background);
441  return diffIm;
442  }
443 
444 /***********************************************************************************************************/
445 //
446 // Explicit instantiations
447 //
448  typedef float PixelT;
449 
450  template class KernelCandidate<PixelT>;
451 
452 }}} // end of namespace lsst::ip::diffim
An include file to include the public header files for lsst::afw::math.
boost::shared_ptr< ImageT const > getImage() const
void setStatus(Status status)
Set the candidate&#39;s status.
Definition: SpatialCell.cc:61
int getInt(const std::string &name) const
return an integer value associated with the given name.
Definition: Policy.h:603
void apply(lsst::afw::image::MaskedImage< PixelT > const &image)
MaskedImagePtr _scienceMaskedImage
Subimage around which you build kernel.
Class stored in SpatialCells for spatial Kernel fitting.
a container for holding hierarchical configuration data in memory.
Definition: Policy.h:169
std::shared_ptr< const Image< PixelT > > ConstPtr
Definition: Image.h:420
float getYCenter() const
Return the object&#39;s row-centre.
Definition: SpatialCell.h:98
afw::image::MaskedImage< PixelT > getDifferenceImage(CandidateSwitch cand)
Calculate associated difference image using internal solutions.
Image Subtraction helper functions.
void build(afw::math::KernelList const &basisList)
Core functionality of KernelCandidate, to build and fill a KernelSolution.
std::shared_ptr< StaticKernelSolution< PixelT > > getKernelSolution(CandidateSwitch cand) const
afw::math::Kernel::Ptr getKernel(CandidateSwitch cand) const
Return results of kernel solution.
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:513
void _buildKernelSolution(afw::math::KernelList const &basisList, std::shared_ptr< Eigen::MatrixXd > hMat)
A class to evaluate image statistics.
Definition: Statistics.h:212
afwImage::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)
Implement fundamental difference imaging step of convolution and subtraction : D = I - (K*T + bg) whe...
Declaration of classes to store the solution for convolution kernels.
VariancePtr getVariance(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s variance.
Definition: MaskedImage.h:896
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
LSST DM logging module built on log4cxx.
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1034
An include file to include the header files for lsst::afw::image.
boost::shared_ptr< Kernel > Ptr
Definition: Kernel.h:138
std::shared_ptr< afw::image::Image< afw::image::VariancePixel > > VariancePtr
estimate sample median
Definition: Statistics.h:70
pex::policy::Policy _policy
Policy.
int getId() const
Return the candidate&#39;s unique ID.
Definition: SpatialCell.h:109
double getKsum(CandidateSwitch cand) const
std::shared_ptr< afw::image::MaskedImage< PixelT > > MaskedImagePtr
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:78
float getXCenter() const
Return the object&#39;s column-centre.
Definition: SpatialCell.h:95
std::shared_ptr< Image< PixelT > > Ptr
Definition: Image.h:419
Class used by SpatialModelCell for spatial Kernel fitting.
Class to calculate difference image statistics.
double _coreFlux
Mean S/N in the science image.
#define LSST_EXCEPT(type,...)
Create an exception with a given type and message and optionally other arguments (dependent on the ty...
Definition: Exception.h:46
double getBackground(CandidateSwitch cand) const
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
Definition: Statistics.cc:1107
boost::shared_ptr< ImageT > getKernelImage(CandidateSwitch cand) const
KernelCandidate(float const xCenter, float const yCenter, MaskedImagePtr const &templateMaskedImage, MaskedImagePtr const &scienceMaskedImage, pex::policy::Policy const &policy)
Constructor.
std::shared_ptr< afw::table::SourceRecord > SourcePtr
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:416
Image Subtraction helper functions.
std::vector< boost::shared_ptr< Kernel > > KernelList
Definition: Kernel.h:539