LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
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 
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
#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:504
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:58
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:1079
float getYCenter() const
Return the object's row-centre.
Definition: SpatialCell.h:93
float getXCenter() const
Return the object's column-centre.
Definition: SpatialCell.h:90
int getId() const
Return the candidate's unique ID.
Definition: SpatialCell.h:104
void setStatus(Status status)
Set the candidate's status.
Definition: SpatialCell.cc:54
A class to evaluate image statistics.
Definition: Statistics.h:215
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1056
Class for storing generic metadata.
Definition: PropertySet.h:67
Class to calculate difference image statistics.
void apply(lsst::afw::image::MaskedImage< PixelT > const &image)
Class stored in SpatialCells for spatial Kernel fitting.
boost::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.
boost::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:354
@ MEDIAN
estimate sample median
Definition: Statistics.h:70
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.