LSST Applications g042eb84c57+730a74494b,g04e9c324dd+8c5ae1fdc5,g134cb467dc+1f1e3e7524,g199a45376c+0ba108daf9,g1fd858c14a+fa7d31856b,g210f2d0738+f66ac109ec,g262e1987ae+83a3acc0e5,g29ae962dfc+d856a2cb1f,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+a1e0c9f713,g47891489e3+0d594cb711,g4d44eb3520+c57ec8f3ed,g4d7b6aa1c5+f66ac109ec,g53246c7159+8c5ae1fdc5,g56a1a4eaf3+fd7ad03fde,g64539dfbff+f66ac109ec,g67b6fd64d1+0d594cb711,g67fd3c3899+f66ac109ec,g6985122a63+0d594cb711,g74acd417e5+3098891321,g786e29fd12+668abc6043,g81db2e9a8d+98e2ab9f28,g87389fa792+8856018cbb,g89139ef638+0d594cb711,g8d7436a09f+80fda9ce03,g8ea07a8fe4+760ca7c3fc,g90f42f885a+033b1d468d,g97be763408+a8a29bda4b,g99822b682c+e3ec3c61f9,g9d5c6a246b+0d5dac0c3d,ga41d0fce20+9243b26dd2,gbf99507273+8c5ae1fdc5,gd7ef33dd92+0d594cb711,gdab6d2f7ff+3098891321,ge410e46f29+0d594cb711,geaed405ab2+c4bbc419c6,gf9a733ac38+8c5ae1fdc5,w.2025.38
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