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
AssessSpatialKernelVisitor.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
12 #include "lsst/afw/math.h"
13 #include "lsst/afw/image.h"
14 #include "lsst/pex/policy/Policy.h"
16 #include "lsst/pex/logging/Trace.h"
17 
21 
22 #define DEBUG_IMAGES 0
23 
24 namespace afwMath = lsst::afw::math;
25 namespace afwImage = lsst::afw::image;
26 namespace pexLogging = lsst::pex::logging;
27 namespace pexPolicy = lsst::pex::policy;
28 namespace pexExcept = lsst::pex::exceptions;
29 
30 namespace lsst {
31 namespace ip {
32 namespace diffim {
33 namespace detail {
54  template<typename PixelT>
58  lsst::pex::policy::Policy const& policy
59  ) :
60  afwMath::CandidateVisitor(),
61  _spatialKernel(spatialKernel),
62  _spatialBackground(spatialBackground),
63  _policy(policy),
64  _imstats(ImageStatistics<PixelT>(_policy)),
65  _nGood(0),
66  _nRejected(0),
67  _nProcessed(0),
68  _useCoreStats(_policy.getBool("useCoreStats")),
69  _coreRadius(_policy.getInt("candidateCoreRadius"))
70  {};
71 
72  template<typename PixelT>
75  ) {
76 
77  KernelCandidate<PixelT> *kCandidate = dynamic_cast<KernelCandidate<PixelT> *>(candidate);
78  if (kCandidate == NULL) {
79  throw LSST_EXCEPT(pexExcept::LogicError,
80  "Failed to cast SpatialCellCandidate to KernelCandidate");
81  }
82  if (!(kCandidate->isInitialized())) {
84  pexLogging::TTrace<3>("lsst.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
85  "Cannot process candidate %d, continuing", kCandidate->getId());
86  return;
87  }
88 
89  pexLogging::TTrace<2>("lsst.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
90  "Processing candidate %d", kCandidate->getId());
91 
92  /*
93  Note - this is a hack until the Kernel API is upgraded by the
94  Davis crew. I need a "local" version of the spatially varying
95  Kernel
96  */
97  afwImage::Image<double> kImage(_spatialKernel->getDimensions());
98  double kSum = _spatialKernel->computeImage(kImage, false,
99  kCandidate->getXCenter(), kCandidate->getYCenter());
100  boost::shared_ptr<afwMath::Kernel>
101  kernelPtr(new afwMath::FixedKernel(kImage));
102  /* </hack> */
103 
104  double background = (*_spatialBackground)(kCandidate->getXCenter(), kCandidate->getYCenter());
105 
106  MaskedImageT diffim = kCandidate->getDifferenceImage(kernelPtr, background);
107 
108  if (DEBUG_IMAGES) {
109  kImage.writeFits(str(boost::format("askv_k%d.fits") % kCandidate->getId()));
110  diffim.writeFits(str(boost::format("askv_d%d.fits") % kCandidate->getId()));
111  }
112 
113  /* Official resids */
114  try {
115  if (_useCoreStats)
116  _imstats.apply(diffim, _coreRadius);
117  else
118  _imstats.apply(diffim);
119  } catch (pexExcept::Exception& e) {
120  pexLogging::TTrace<3>("lsst.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
121  "Unable to calculate imstats for Candidate %d", kCandidate->getId());
123  return;
124  }
125 
126  _nProcessed += 1;
127 
128  pexLogging::TTrace<5>("lsst.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
129  "Chi2 = %.3f", _imstats.getVariance());
130  pexLogging::TTrace<5>("lsst.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
131  "X = %.2f Y = %.2f",
132  kCandidate->getXCenter(),
133  kCandidate->getYCenter());
134  pexLogging::TTrace<5>("lsst.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
135  "Kernel Sum = %.3f", kSum);
136  pexLogging::TTrace<5>("lsst.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
137  "Background = %.3f", background);
138  pexLogging::TTrace<3>("lsst.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
139  "Candidate %d resids = %.3f +/- %.3f sigma (%d pix)",
140  kCandidate->getId(),
141  _imstats.getMean(),
142  _imstats.getRms(),
143  _imstats.getNpix());
144 
145  bool meanIsNan = std::isnan(_imstats.getMean());
146  bool rmsIsNan = std::isnan(_imstats.getRms());
147  if (meanIsNan || rmsIsNan) {
149  pexLogging::TTrace<4>("lsst.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
150  "Rejecting candidate %d, encountered NaN",
151  kCandidate->getId());
152  _nRejected += 1;
153  return;
154  }
155 
156  if (_policy.getBool("spatialKernelClipping")) {
157  if (fabs(_imstats.getMean()) > _policy.getDouble("candidateResidualMeanMax")) {
159  pexLogging::TTrace<4>("lsst.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
160  "Rejecting candidate %d; bad mean residual : |%.3f| > %.3f",
161  kCandidate->getId(),
162  _imstats.getMean(),
163  _policy.getDouble("candidateResidualMeanMax"));
164  _nRejected += 1;
165  }
166  else if (_imstats.getRms() > _policy.getDouble("candidateResidualStdMax")) {
168  pexLogging::TTrace<4>("lsst.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
169  "Rejecting candidate %d; bad residual rms : %.3f > %.3f",
170  kCandidate->getId(),
171  _imstats.getRms(),
172  _policy.getDouble("candidateResidualStdMax"));
173  _nRejected += 1;
174  }
175  else {
177  pexLogging::TTrace<4>("lsst.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
178  "Spatial kernel OK");
179  _nGood += 1;
180  }
181  }
182  else {
184  pexLogging::TTrace<6>("lsst.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
185  "Sigma clipping not enabled");
186  _nGood += 1;
187  }
188 
189  /* Core resids for debugging */
190  if (!(_useCoreStats)) {
191  try {
192  _imstats.apply(diffim, _coreRadius);
193  } catch (pexExcept::Exception& e) {
194  pexLogging::TTrace<3>("lsst.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
195  "Unable to calculate core imstats for Candidate %d",
196  kCandidate->getId());
198  return;
199  }
200  pexLogging::TTrace<4>("lsst.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
201  "Candidate %d core resids = %.3f +/- %.3f sigma (%d pix)",
202  kCandidate->getId(),
203  _imstats.getMean(),
204  _imstats.getRms(),
205  _imstats.getNpix());
206  }
207  }
208 
209  typedef float PixelT;
210  template class AssessSpatialKernelVisitor<PixelT>;
211 
212 }}}} // end of namespace lsst::ip::diffim::detail
An include file to include the public header files for lsst::afw::math.
boost::shared_ptr< lsst::afw::math::Function2< double > > SpatialFunctionPtr
Definition: Kernel.h:143
Asseses the quality of a candidate given a spatial kernel and background model.
float getYCenter() const
Return the object&#39;s row-centre.
Definition: SpatialCell.h:98
Class stored in SpatialCells for spatial Kernel fitting.
AssessSpatialKernelVisitor(lsst::afw::math::LinearCombinationKernel::Ptr spatialKernel, lsst::afw::math::Kernel::SpatialFunctionPtr spatialBackground, lsst::pex::policy::Policy const &policy)
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.
void setStatus(Status status)
Set the candidate&#39;s status.
Definition: SpatialCell.cc:61
int isnan(T t)
Definition: ieee.h:110
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
An include file to include the header files for lsst::afw::image.
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.
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
void processCandidate(lsst::afw::math::SpatialCellCandidate *candidate)
Declaration of AssessSpatialKernelVisitor.
A class to represent a 2-dimensional array of pixels.
Definition: PSF.h:43
Image Subtraction helper functions.
boost::shared_ptr< LinearCombinationKernel > Ptr
Definition: Kernel.h:818
A kernel created from an Image.
Definition: Kernel.h:551
#define DEBUG_IMAGES
float getXCenter() const
Return the object&#39;s column-centre.
Definition: SpatialCell.h:95