LSST Applications g04e9c324dd+8c5ae1fdc5,g134cb467dc+b203dec576,g18429d2f64+358861cd2c,g199a45376c+0ba108daf9,g1fd858c14a+dd066899e3,g262e1987ae+ebfced1d55,g29ae962dfc+72fd90588e,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+b668f15bc5,g4595892280+3897dae354,g47891489e3+abcf9c3559,g4d44eb3520+fb4ddce128,g53246c7159+8c5ae1fdc5,g67b6fd64d1+abcf9c3559,g67fd3c3899+1f72b5a9f7,g74acd417e5+cb6b47f07b,g786e29fd12+668abc6043,g87389fa792+8856018cbb,g89139ef638+abcf9c3559,g8d7436a09f+bcf525d20c,g8ea07a8fe4+9f5ccc88ac,g90f42f885a+6054cc57f1,g97be763408+06f794da49,g9dd6db0277+1f72b5a9f7,ga681d05dcb+7e36ad54cd,gabf8522325+735880ea63,gac2eed3f23+abcf9c3559,gb89ab40317+abcf9c3559,gbf99507273+8c5ae1fdc5,gd8ff7fe66e+1f72b5a9f7,gdab6d2f7ff+cb6b47f07b,gdc713202bf+1f72b5a9f7,gdfd2d52018+8225f2b331,ge365c994fd+375fc21c71,ge410e46f29+abcf9c3559,geaed405ab2+562b3308c0,gf9a733ac38+8c5ae1fdc5,w.2025.35
LSST Data Management Base Package
Loading...
Searching...
No Matches
AssessSpatialKernelVisitor.cc
Go to the documentation of this file.
1// -*- lsst-c++ -*-
11
12#include "lsst/afw/math.h"
13#include "lsst/afw/image.h"
14#include "lsst/log/Log.h"
17
21
22#define DEBUG_IMAGES 0
23
24namespace afwMath = lsst::afw::math;
25namespace afwImage = lsst::afw::image;
26namespace dafBase = lsst::daf::base;
27namespace pexExcept = lsst::pex::exceptions;
28
29namespace lsst {
30namespace ip {
31namespace diffim {
32namespace detail {
53 template<typename PixelT>
58 ) :
60 _spatialKernel(spatialKernel),
61 _spatialBackground(spatialBackground),
62 _ps(ps.deepCopy()),
63 _imstats(ImageStatistics<PixelT>(ps)),
64 _nGood(0),
65 _nRejected(0),
66 _nProcessed(0),
67 _useCoreStats(ps.getAsBool("useCoreStats")),
68 _coreRadius(ps.getAsInt("candidateCoreRadius"))
69 {};
70
71 template<typename PixelT>
74 ) {
75
76 KernelCandidate<PixelT> *kCandidate = dynamic_cast<KernelCandidate<PixelT> *>(candidate);
77 if (kCandidate == NULL) {
78 throw LSST_EXCEPT(pexExcept::LogicError,
79 "Failed to cast SpatialCellCandidate to KernelCandidate");
80 }
81 if (!(kCandidate->isInitialized())) {
83 LOGL_DEBUG("TRACE2.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
84 "Cannot process candidate %d, continuing", kCandidate->getId());
85 return;
86 }
87
88 LOGL_DEBUG("TRACE1.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
89 "Processing candidate %d", kCandidate->getId());
90
91 /*
92 Note - this is a hack until the Kernel API is upgraded by the
93 Davis crew. I need a "local" version of the spatially varying
94 Kernel
95 */
96 afwImage::Image<double> kImage(_spatialKernel->getDimensions());
97 double kSum = _spatialKernel->computeImage(kImage, false,
98 kCandidate->getXCenter(), kCandidate->getYCenter());
100 kernelPtr(new afwMath::FixedKernel(kImage));
101 /* </hack> */
102
103 double background = (*_spatialBackground)(kCandidate->getXCenter(), kCandidate->getYCenter());
104
105 MaskedImageT diffim = kCandidate->getDifferenceImage(kernelPtr, background);
106
107 if (DEBUG_IMAGES) {
108 kImage.writeFits(str(boost::format("askv_k%d.fits") % kCandidate->getId()));
109 diffim.writeFits(str(boost::format("askv_d%d.fits") % kCandidate->getId()));
110 }
111
112 /* Official resids */
113 try {
114 if (_useCoreStats)
115 _imstats.apply(diffim, _coreRadius);
116 else
117 _imstats.apply(diffim);
118 } catch (pexExcept::Exception& e) {
119 LOGL_DEBUG("TRACE2.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
120 "Unable to calculate imstats for Candidate %d", kCandidate->getId());
122 return;
123 }
124
125 _nProcessed += 1;
126
127 LOGL_DEBUG("TRACE4.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
128 "Chi2 = %.3f", _imstats.getVariance());
129 LOGL_DEBUG("TRACE4.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
130 "X = %.2f Y = %.2f",
131 kCandidate->getXCenter(),
132 kCandidate->getYCenter());
133 LOGL_DEBUG("TRACE4.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
134 "Kernel Sum = %.3f", kSum);
135 LOGL_DEBUG("TRACE4.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
136 "Background = %.3f", background);
137 LOGL_DEBUG("TRACE2.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
138 "Candidate %d resids = %.3f +/- %.3f sigma (%d pix)",
139 kCandidate->getId(),
140 _imstats.getMean(),
141 _imstats.getRms(),
142 _imstats.getNpix());
143
144 bool meanIsNan = std::isnan(_imstats.getMean());
145 bool rmsIsNan = std::isnan(_imstats.getRms());
146 if (meanIsNan || rmsIsNan) {
148 LOGL_DEBUG("TRACE3.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
149 "Rejecting candidate %d, encountered NaN",
150 kCandidate->getId());
151 _nRejected += 1;
152 return;
153 }
154
155 if (_ps->getAsBool("spatialKernelClipping")) {
156 if (fabs(_imstats.getMean()) > _ps->getAsDouble("candidateResidualMeanMax")) {
158 LOGL_DEBUG("TRACE3.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
159 "Rejecting candidate %d; bad mean residual : |%.3f| > %.3f",
160 kCandidate->getId(),
161 _imstats.getMean(),
162 _ps->getAsDouble("candidateResidualMeanMax"));
163 _nRejected += 1;
164 }
165 else if (_imstats.getRms() > _ps->getAsDouble("candidateResidualStdMax")) {
167 LOGL_DEBUG("TRACE3.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
168 "Rejecting candidate %d; bad residual rms : %.3f > %.3f",
169 kCandidate->getId(),
170 _imstats.getRms(),
171 _ps->getAsDouble("candidateResidualStdMax"));
172 _nRejected += 1;
173 }
174 else {
176 LOGL_DEBUG("TRACE3.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
177 "Spatial kernel OK");
178 _nGood += 1;
179 }
180 }
181 else {
183 LOGL_DEBUG("TRACE5.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
184 "Sigma clipping not enabled");
185 _nGood += 1;
186 }
187
188 /* Core resids for debugging */
189 if (!(_useCoreStats)) {
190 try {
191 _imstats.apply(diffim, _coreRadius);
192 } catch (pexExcept::Exception& e) {
193 LOGL_DEBUG("TRACE2.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
194 "Unable to calculate core imstats for Candidate %d",
195 kCandidate->getId());
197 return;
198 }
199 LOGL_DEBUG("TRACE3.ip.diffim.AssessSpatialKernelVisitor.processCandidate",
200 "Candidate %d core resids = %.3f +/- %.3f sigma (%d pix)",
201 kCandidate->getId(),
202 _imstats.getMean(),
203 _imstats.getRms(),
204 _imstats.getNpix());
205 }
206 }
207
208 typedef float PixelT;
210
211}}}} // end of namespace lsst::ip::diffim::detail
Declaration of AssessSpatialKernelVisitor.
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
Image Subtraction helper functions.
Class used by SpatialModelCell for spatial Kernel fitting.
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:515
A class to represent a 2-dimensional array of pixels.
Definition Image.h:51
void writeFits(std::string const &fileName, daf::base::PropertySet const *metadata=nullptr, std::string const &mode="w") const
Write an image to a regular FITS file.
A kernel created from an Image.
Definition Kernel.h:471
std::shared_ptr< lsst::afw::math::Function2< double > > SpatialFunctionPtr
Definition Kernel.h:113
Base class for candidate objects in a SpatialCell.
Definition SpatialCell.h:70
float getYCenter() const
Return the object's row-centre.
Definition SpatialCell.h:91
float getXCenter() const
Return the object's column-centre.
Definition SpatialCell.h:88
int getId() const
Return the candidate's unique ID.
void setStatus(Status status)
Set the candidate's status.
Class for storing generic metadata.
Definition PropertySet.h:67
Class to calculate difference image statistics.
Class stored in SpatialCells for spatial Kernel fitting.
afw::image::MaskedImage< PixelT > getDifferenceImage(CandidateSwitch cand)
Calculate associated difference image using internal solutions.
Asseses the quality of a candidate given a spatial kernel and background model.
AssessSpatialKernelVisitor(std::shared_ptr< lsst::afw::math::LinearCombinationKernel > spatialKernel, lsst::afw::math::Kernel::SpatialFunctionPtr spatialBackground, lsst::daf::base::PropertySet const &ps)
void processCandidate(lsst::afw::math::SpatialCellCandidate *candidate)
T isnan(T... args)
#define DEBUG_IMAGES