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
BuildSingleKernelVisitor.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
12 #include "boost/shared_ptr.hpp"
13 #include "Eigen/Core"
14 
15 #include "lsst/afw/math.h"
16 #include "lsst/afw/image.h"
17 
19 #include "lsst/pex/policy/Policy.h"
20 #include "lsst/pex/logging/Trace.h"
21 
25 
26 #define DEBUG_MATRIX 0
27 
28 namespace afwMath = lsst::afw::math;
29 namespace afwImage = lsst::afw::image;
30 namespace pexLogging = lsst::pex::logging;
31 namespace pexPolicy = lsst::pex::policy;
32 namespace pexExcept = lsst::pex::exceptions;
33 namespace ipDiffim = lsst::ip::diffim;
34 
35 namespace lsst {
36 namespace ip {
37 namespace diffim {
38 namespace detail {
39 
89  template<typename PixelT>
91  lsst::afw::math::KernelList const& basisList,
92  lsst::pex::policy::Policy const& policy
94  ) :
95  afwMath::CandidateVisitor(),
96  _basisList(basisList),
97  _policy(policy),
98  _hMat(),
99  _imstats(ImageStatistics<PixelT>(_policy)),
100  _skipBuilt(true),
101  _nRejected(0),
102  _nProcessed(0),
103  _useRegularization(false),
104  _useCoreStats(_policy.getBool("useCoreStats")),
105  _coreRadius(_policy.getInt("candidateCoreRadius"))
106  {};
107 
108  template<typename PixelT>
110  lsst::afw::math::KernelList const& basisList,
111  lsst::pex::policy::Policy const& policy,
113  boost::shared_ptr<Eigen::MatrixXd> hMat
114  ) :
115  afwMath::CandidateVisitor(),
116  _basisList(basisList),
117  _policy(policy),
118  _hMat(hMat),
119  _imstats(ImageStatistics<PixelT>(_policy)),
120  _skipBuilt(true),
121  _nRejected(0),
122  _nProcessed(0),
123  _useRegularization(true),
124  _useCoreStats(_policy.getBool("useCoreStats")),
125  _coreRadius(_policy.getInt("candidateCoreRadius"))
126  {};
127 
128 
129  template<typename PixelT>
132  ) {
133 
134  ipDiffim::KernelCandidate<PixelT> *kCandidate =
135  dynamic_cast<ipDiffim::KernelCandidate<PixelT> *>(candidate);
136  if (kCandidate == NULL) {
137  pexLogging::TTrace<3>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
138  "Failed to cast SpatialCellCandidate to KernelCandidate %d",
139  kCandidate->getId());
140  throw LSST_EXCEPT(pexExcept::LogicError,
141  "Failed to cast SpatialCellCandidate to KernelCandidate");
142  }
143 
144  if (_skipBuilt and kCandidate->isInitialized()) {
145  return;
146  }
147 
148  pexLogging::TTrace<2>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
149  "Processing candidate %d", kCandidate->getId());
150  pexLogging::TTrace<5>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
151  "X = %.2f Y = %.2f",
152  kCandidate->getXCenter(),
153  kCandidate->getYCenter());
154 
155  /* Build its kernel here */
156  try {
157  if (_useRegularization)
158  kCandidate->build(_basisList, _hMat);
159  else
160  kCandidate->build(_basisList);
161 
162  } catch (pexExcept::Exception &e) {
164  pexLogging::TTrace<4>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
165  "Unable to process candidate %d; exception caught (%s)",
166  kCandidate->getId(),
167  e.what());
168  _nRejected += 1;
169  return;
170  }
171 
172  if (kCandidate->getStatus() == afwMath::SpatialCellCandidate::BAD) {
173  pexLogging::TTrace<4>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
174  "Candidate %d Returned BAD upon build, exiting",
175  kCandidate->getId());
176  _nRejected += 1;
177  return;
178  }
179 
180 
181  /*
182  * Make diffim and set chi2 from result. Note that you need to use the
183  * most recent kernel
184  */
186  try {
187  if (_useCoreStats)
188  _imstats.apply(diffim, _coreRadius);
189  else
190  _imstats.apply(diffim);
191  } catch (pexExcept::Exception& e) {
192  pexLogging::TTrace<3>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
193  "Unable to calculate imstats for Candidate %d", kCandidate->getId());
195  return;
196  }
197  _nProcessed += 1;
198 
199  kCandidate->setChi2(_imstats.getVariance());
200 
201  /* When using a Pca basis, we don't reset the kernel or background,
202  so we need to evaluate these locally for the Trace */
203  double kSum = kCandidate->getKsum(ipDiffim::KernelCandidate<PixelT>::RECENT);
204  double background = kCandidate->getBackground(ipDiffim::KernelCandidate<PixelT>::RECENT);
205 
206  pexLogging::TTrace<5>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
207  "Chi2 = %.3f", kCandidate->getChi2());
208  pexLogging::TTrace<5>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
209  "Kernel Sum = %.3f", kSum);
210  pexLogging::TTrace<5>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
211  "Background = %.3f", background);
212  pexLogging::TTrace<3>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
213  "Candidate %d resids = %.3f +/- %.3f sigma (%d pix)",
214  kCandidate->getId(),
215  _imstats.getMean(),
216  _imstats.getRms(),
217  _imstats.getNpix());
218 
219  bool meanIsNan = std::isnan(_imstats.getMean());
220  bool rmsIsNan = std::isnan(_imstats.getRms());
221  if (meanIsNan || rmsIsNan) {
223  pexLogging::TTrace<4>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
224  "Rejecting candidate %d, encountered NaN",
225  kCandidate->getId());
226  _nRejected += 1;
227  return;
228  }
229 
230  if (_policy.getBool("singleKernelClipping")) {
231  if (fabs(_imstats.getMean()) > _policy.getDouble("candidateResidualMeanMax")) {
233  pexLogging::TTrace<4>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
234  "Rejecting candidate %d; bad mean residual : |%.3f| > %.3f",
235  kCandidate->getId(),
236  _imstats.getMean(),
237  _policy.getDouble("candidateResidualMeanMax"));
238  _nRejected += 1;
239  }
240  else if (_imstats.getRms() > _policy.getDouble("candidateResidualStdMax")) {
242  pexLogging::TTrace<4>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
243  "Rejecting candidate %d; bad residual rms : %.3f > %.3f",
244  kCandidate->getId(),
245  _imstats.getRms(),
246  _policy.getDouble("candidateResidualStdMax"));
247  _nRejected += 1;
248  }
249  else {
251  pexLogging::TTrace<4>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
252  "Source kernel OK");
253  }
254  }
255  else {
257  pexLogging::TTrace<6>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
258  "Sigma clipping not enabled");
259  }
260 
261  /* Core resids for debugging */
262  if (!(_useCoreStats)) {
263  try {
264  _imstats.apply(diffim, _coreRadius);
265  } catch (pexExcept::Exception& e) {
266  pexLogging::TTrace<3>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
267  "Unable to calculate core imstats for Candidate %d",
268  kCandidate->getId());
270  return;
271  }
272  pexLogging::TTrace<4>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
273  "Candidate %d core resids = %.3f +/- %.3f sigma (%d pix)",
274  kCandidate->getId(),
275  _imstats.getMean(),
276  _imstats.getRms(),
277  _imstats.getNpix());
278  }
279 
280  }
281 
282  typedef float PixelT;
283 
284  template class BuildSingleKernelVisitor<PixelT>;
285 
286  template boost::shared_ptr<BuildSingleKernelVisitor<PixelT> >
289 
290  template boost::shared_ptr<BuildSingleKernelVisitor<PixelT> >
293  boost::shared_ptr<Eigen::MatrixXd>);
294 
295 }}}} // end of namespace lsst::ip::diffim::detail
An include file to include the public header files for lsst::afw::math.
float getYCenter() const
Return the object&#39;s row-centre.
Definition: SpatialCell.h:98
Class stored in SpatialCells for spatial Kernel fitting.
void processCandidate(lsst::afw::math::SpatialCellCandidate *candidate)
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
Builds the convolution kernel for a given candidate.
afw::image::MaskedImage< PixelT > getDifferenceImage(CandidateSwitch cand)
Calculate associated difference image using internal solutions.
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
template boost::shared_ptr< BuildSingleKernelVisitor< PixelT > > makeBuildSingleKernelVisitor< PixelT >(lsst::afw::math::KernelList const &, lsst::pex::policy::Policy const &)
int isnan(T t)
Definition: ieee.h:110
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
BuildSingleKernelVisitor(lsst::afw::math::KernelList const &basisList, lsst::pex::policy::Policy const &policy)
An include file to include the header files for lsst::afw::image.
Status getStatus() const
Return the candidate&#39;s status.
Definition: SpatialCell.h:111
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.
void setChi2(double chi2)
Set the candidate&#39;s chi^2.
Definition: SpatialCell.h:165
virtual char const * what(void) const
Definition: Exception.cc:112
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
double getChi2() const
Return the candidate&#39;s chi^2.
Definition: SpatialCell.h:163
double getBackground(CandidateSwitch cand) const
std::vector< boost::shared_ptr< Kernel > > KernelList
Definition: Kernel.h:542
Implementation of BuildSingleKernelVisitor.
Image Subtraction helper functions.
float getXCenter() const
Return the object&#39;s column-centre.
Definition: SpatialCell.h:95