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