LSSTApplications  17.0+11,17.0+34,17.0+56,17.0+57,17.0+59,17.0+7,17.0-1-g377950a+33,17.0.1-1-g114240f+2,17.0.1-1-g4d4fbc4+28,17.0.1-1-g55520dc+49,17.0.1-1-g5f4ed7e+52,17.0.1-1-g6dd7d69+17,17.0.1-1-g8de6c91+11,17.0.1-1-gb9095d2+7,17.0.1-1-ge9fec5e+5,17.0.1-1-gf4e0155+55,17.0.1-1-gfc65f5f+50,17.0.1-1-gfc6fb1f+20,17.0.1-10-g87f9f3f+1,17.0.1-11-ge9de802+16,17.0.1-16-ga14f7d5c+4,17.0.1-17-gc79d625+1,17.0.1-17-gdae4c4a+8,17.0.1-2-g26618f5+29,17.0.1-2-g54f2ebc+9,17.0.1-2-gf403422+1,17.0.1-20-g2ca2f74+6,17.0.1-23-gf3eadeb7+1,17.0.1-3-g7e86b59+39,17.0.1-3-gb5ca14a,17.0.1-3-gd08d533+40,17.0.1-30-g596af8797,17.0.1-4-g59d126d+4,17.0.1-4-gc69c472+5,17.0.1-6-g5afd9b9+4,17.0.1-7-g35889ee+1,17.0.1-7-gc7c8782+18,17.0.1-9-gc4bbfb2+3,w.2019.22
LSSTDataManagementBasePackage
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 
20 namespace afwMath = lsst::afw::math;
21 namespace afwImage = lsst::afw::image;
22 namespace pexExcept = lsst::pex::exceptions;
23 
24 namespace lsst {
25 namespace ip {
26 namespace diffim {
27 namespace 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);
81  new afwMath::FixedKernel(img)
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  PTR(ImageT) kImage = kCandidate->getKernelSolution(
102  KernelCandidate<PixelT>::ORIG)->makeKernelImage();
103  *kImage /= kCandidate->getKernelSolution(
104  KernelCandidate<PixelT>::ORIG)->getKsum();
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  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
uint64_t * ptr
Definition: RangeSet.cc:88
estimate sample minimum
Definition: Statistics.h:76
Class stored in SpatialCells for spatial Kernel fitting.
KernelPcaVisitor(std::shared_ptr< KernelPca< ImageT > > imagePca)
Definition: KernelPca.cc:56
std::shared_ptr< StaticKernelSolution< PixelT > > getKernelSolution(CandidateSwitch cand) const
estimate sample maximum
Definition: Statistics.h:77
T end(T... args)
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
int getId() const
Return the candidate&#39;s unique ID.
Definition: SpatialCell.h:104
A class to run a PCA on all candidate kernels (represented as Images).
Definition: KernelPca.h:40
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:489
int min
A class to evaluate image statistics.
Definition: Statistics.h:215
void processCandidate(lsst::afw::math::SpatialCellCandidate *candidate)
Definition: KernelPca.cc:89
LSST DM logging module built on log4cxx.
lsst::afw::math::KernelList getEigenKernels()
Definition: KernelPca.cc:65
T push_back(T... args)
A base class for image defects.
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:354
Overrides the analyze method of base class afwImage::ImagePca.
Definition: KernelPca.h:24
int max
Declaration of KernelPca and KernelPcaVisitor.
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
Class used by SpatialModelCell for spatial Kernel fitting.
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1056
T size(T... args)
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
T begin(T... args)
afw::table::Key< afw::table::Array< ImagePixelT > > image
#define PTR(...)
Definition: base.h:41
Base class for candidate objects in a SpatialCell.
Definition: SpatialCell.h:70
lsst::afw::image::Image< lsst::afw::math::Kernel::Pixel > ImageT
Definition: KernelPca.h:42
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:59
int end
virtual void analyze()
Generate eigenimages that are normalised.
Definition: KernelPca.cc:163
A kernel created from an Image.
Definition: Kernel.h:509