LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
KernelPca.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/log/Log.h"
16
19
20namespace afwMath = lsst::afw::math;
21namespace afwImage = lsst::afw::image;
23
24namespace lsst {
25namespace ip {
26namespace diffim {
27namespace detail {
28
55 template<typename PixelT>
58 ) :
59 afwMath::CandidateVisitor(),
60 _imagePca(imagePca),
61 _mean()
62 {};
63
64 template<typename PixelT>
66 afwMath::KernelList kernelList;
67
68 std::vector<std::shared_ptr<ImageT>> eigenImages = _imagePca->getEigenImages();
69 int ncomp = eigenImages.size();
70
71 if (_mean) {
75 (*_mean, true))));
76 }
77 for (int i = 0; i < ncomp; i++) {
79 afwImage::Image<afwMath::Kernel::Pixel>(*eigenImages[i], true);
82 ));
83 }
84
85 return kernelList;
86 }
87
88 template<typename PixelT>
90
91 KernelCandidate<PixelT> *kCandidate = dynamic_cast<KernelCandidate<PixelT> *>(candidate);
92 if (kCandidate == NULL) {
94 "Failed to cast SpatialCellCandidate to KernelCandidate");
95 }
96 LOGL_DEBUG("TRACE5.ip.diffim.SetPcaImageVisitor.processCandidate",
97 "Processing candidate %d", kCandidate->getId());
98
99 try {
100 /* Normalize to unit sum */
101 std::shared_ptr<ImageT> kImage = kCandidate->getKernelSolution(
102 KernelCandidate<PixelT>::ORIG)->makeKernelImage();
103 *kImage /= kCandidate->getKernelSolution(
105 /* Tell imagePca they have the same weighting in the Pca */
106 _imagePca->addImage(kImage, 1.0);
107 } catch(pexExcept::Exception &e) {
108 return;
109 }
110 }
111
112 template<typename PixelT>
114 /*
115 If we don't subtract off the mean before we do the Pca, the
116 subsequent terms carry less of the power than if you do subtract
117 off the mean. Explicit example:
118
119 With mean subtraction:
120 DEBUG: Eigenvalue 0 : 0.010953 (0.373870 %)
121 DEBUG: Eigenvalue 1 : 0.007927 (0.270604 %)
122 DEBUG: Eigenvalue 2 : 0.001393 (0.047542 %)
123 DEBUG: Eigenvalue 3 : 0.001092 (0.037261 %)
124 DEBUG: Eigenvalue 4 : 0.000829 (0.028283 %)
125
126 Without mean subtraction:
127 DEBUG: Eigenvalue 0 : 0.168627 (0.876046 %)
128 DEBUG: Eigenvalue 1 : 0.007935 (0.041223 %)
129 DEBUG: Eigenvalue 2 : 0.006049 (0.031424 %)
130 DEBUG: Eigenvalue 3 : 0.001188 (0.006173 %)
131 DEBUG: Eigenvalue 4 : 0.001050 (0.005452 %)
132
133 After the first term above, which basically represents the mean,
134 the remaining terms carry less of the power than if you do
135 subtract off the mean. (0.041223/(1-0.876046) < 0.373870).
136 */
137 LOGL_DEBUG("TRACE5.ip.diffim.KernelPcaVisitor.subtractMean",
138 "Subtracting mean feature before Pca");
139
140 _mean = _imagePca->getMean();
141 KernelPca<ImageT>::ImageList imageList = _imagePca->getImageList();
142 for (typename KernelPca<ImageT>::ImageList::const_iterator ptr = imageList.begin(),
143 end = imageList.end(); ptr != end; ++ptr) {
144 **ptr -= *_mean;
145 }
146 }
147
162 template <typename ImageT>
164 {
165 Super::analyze();
166
167 typename Super::ImageList const &eImageList = this->getEigenImages();
168 typename Super::ImageList::const_iterator iter = eImageList.begin(), end = eImageList.end();
169 for (size_t i = 0; iter != end; ++i, ++iter) {
170 std::shared_ptr<ImageT> eImage = *iter;
171
172 /*
173 * Normalise eigenImages to have a maximum of 1.0. For n > 0 they
174 * (should) have mean == 0, so we can't use that to normalize
175 */
177 double const min = stats.getValue(afwMath::MIN);
178 double const max = stats.getValue(afwMath::MAX);
179
180 double const extreme = (fabs(min) > max) ? min :max;
181 if (extreme != 0.0) {
182 *eImage /= extreme;
183 }
184 }
185 }
186
187
188 typedef float PixelT;
189 template class KernelPcaVisitor<PixelT>;
191
192}}}} // end of namespace lsst::ip::diffim::detail
int min
int end
int max
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
Class used by SpatialModelCell for spatial Kernel fitting.
Declaration of KernelPca and KernelPcaVisitor.
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
uint64_t * ptr
Definition RangeSet.cc:95
T begin(T... args)
A class to represent a 2-dimensional array of pixels.
Definition Image.h:51
A kernel created from an Image.
Definition Kernel.h:471
Base class for candidate objects in a SpatialCell.
Definition SpatialCell.h:70
int getId() const
Return the candidate's unique ID.
A class to evaluate image statistics.
Definition Statistics.h:222
Class stored in SpatialCells for spatial Kernel fitting.
std::shared_ptr< StaticKernelSolution< PixelT > > getKernelSolution(CandidateSwitch cand) const
Overrides the analyze method of base class afwImage::ImagePca.
Definition KernelPca.h:24
virtual void analyze()
Generate eigenimages that are normalised.
Definition KernelPca.cc:163
A class to run a PCA on all candidate kernels (represented as Images).
Definition KernelPca.h:40
void processCandidate(lsst::afw::math::SpatialCellCandidate *candidate)
Definition KernelPca.cc:89
lsst::afw::math::KernelList getEigenKernels()
Definition KernelPca.cc:65
KernelPcaVisitor(std::shared_ptr< KernelPca< ImageT > > imagePca)
Definition KernelPca.cc:56
Provides consistent interface for LSST exceptions.
Definition Exception.h:107
Reports errors in the logical structure of the program.
Definition Runtime.h:46
T end(T... args)
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:361
@ MIN
estimate sample minimum
Definition Statistics.h:66
@ MAX
estimate sample maximum
Definition Statistics.h:67
T push_back(T... args)
T size(T... args)