LSSTApplications  18.1.0
LSSTDataManagementBasePackage
ImagePca.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
25 /*
26  * Utilities to support PCA analysis of a set of images
27  */
28 #include <algorithm>
29 #include <cmath>
30 #include <cstdint>
31 #include <memory>
32 
33 #include "Eigen/Core"
34 #include "Eigen/SVD"
35 #include "Eigen/Eigenvalues"
36 
39 
40 namespace afwMath = lsst::afw::math;
41 
42 namespace lsst {
43 namespace afw {
44 namespace image {
45 
46 template <typename ImageT>
47 ImagePca<ImageT>::ImagePca(bool constantWeight)
48  : _imageList(),
49  _fluxList(),
50  _dimensions(0, 0),
51  _constantWeight(constantWeight),
52  _eigenValues(std::vector<double>()),
53  _eigenImages(ImageList()) {}
54 
55 template <typename ImageT>
56 ImagePca<ImageT>::ImagePca(ImagePca const&) = default;
57 template <typename ImageT>
59 template <typename ImageT>
61 template <typename ImageT>
63 
64 template <typename ImageT>
66  if (_imageList.empty()) {
67  _dimensions = img->getDimensions();
68  } else {
69  if (getDimensions() != img->getDimensions()) {
71  (boost::format("Dimension mismatch: saw %dx%d; expected %dx%d") %
72  img->getWidth() % img->getHeight() % _dimensions.getX() % _dimensions.getY())
73  .str());
74  }
75  }
76 
77  if (flux == 0.0) {
78  throw LSST_EXCEPT(lsst::pex::exceptions::OutOfRangeError, "Flux may not be zero");
79  }
80 
81  _imageList.push_back(img);
82  _fluxList.push_back(flux);
83 }
84 
85 template <typename ImageT>
87  return _imageList;
88 }
89 
90 template <typename ImageT>
92  if (_imageList.empty()) {
93  throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "You haven't provided any images");
94  }
95 
96  std::shared_ptr<ImageT> mean(new ImageT(getDimensions()));
97  *mean = static_cast<typename ImageT::Pixel>(0);
98 
99  for (typename ImageList::const_iterator ptr = _imageList.begin(), end = _imageList.end(); ptr != end;
100  ++ptr) {
101  *mean += **ptr;
102  }
103  *mean /= _imageList.size();
104 
105  return mean;
106 }
107 
108 namespace {
109 /*
110  * Analyze the images in an ImagePca, calculating the PCA decomposition (== Karhunen-Lo\`eve basis)
111  *
112  * The notation is that in chapter 7 of Gyula Szokoly's thesis at JHU
113  */
114 template <typename T>
115 struct SortEvalueDecreasing
116  : public std::binary_function<std::pair<T, int> const&, std::pair<T, int> const&, bool> {
117  bool operator()(std::pair<T, int> const& a, std::pair<T, int> const& b) {
118  return a.first > b.first; // N.b. sort on greater
119  }
120 };
121 } // namespace
122 
123 template <typename ImageT>
125  int const nImage = _imageList.size();
126 
127  if (nImage == 0) {
128  throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "No images provided for PCA analysis");
129  }
130  /*
131  * Eigen doesn't like 1x1 matrices, but we don't really need it to handle a single matrix...
132  */
133  if (nImage == 1) {
134  _eigenImages.clear();
135  _eigenImages.push_back(std::shared_ptr<ImageT>(new ImageT(*_imageList[0], true)));
136 
137  _eigenValues.clear();
138  _eigenValues.push_back(1.0);
139 
140  return;
141  }
142  /*
143  * Find the eigenvectors/values of the scalar product matrix, R' (Eq. 7.4)
144  */
145  Eigen::MatrixXd R(nImage, nImage); // residuals' inner products
146 
147  double flux_bar = 0; // mean of flux for all regions
148  for (int i = 0; i != nImage; ++i) {
149  typename GetImage<ImageT>::type const& im_i = *GetImage<ImageT>::getImage(_imageList[i]);
150  double const flux_i = getFlux(i);
151  flux_bar += flux_i;
152 
153  for (int j = i; j != nImage; ++j) {
154  typename GetImage<ImageT>::type const& im_j = *GetImage<ImageT>::getImage(_imageList[j]);
155  double const flux_j = getFlux(j);
156 
157  double dot = innerProduct(im_i, im_j);
158  if (_constantWeight) {
159  dot /= flux_i * flux_j;
160  }
161  R(i, j) = R(j, i) = dot / nImage;
162  }
163  }
164  flux_bar /= nImage;
165  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(R);
166  Eigen::MatrixXd const& Q = eVecValues.eigenvectors();
167  Eigen::VectorXd const& lambda = eVecValues.eigenvalues();
168  //
169  // We need to sort the eigenValues, and remember the permutation we applied to the eigenImages
170  // We'll use the vector lambdaAndIndex to achieve this
171  //
172  std::vector<std::pair<double, int> > lambdaAndIndex; // pairs (eValue, index)
173  lambdaAndIndex.reserve(nImage);
174 
175  for (int i = 0; i != nImage; ++i) {
176  lambdaAndIndex.push_back(std::make_pair(lambda(i), i));
177  }
178  std::sort(lambdaAndIndex.begin(), lambdaAndIndex.end(), SortEvalueDecreasing<double>());
179  //
180  // Save the (sorted) eigen values
181  //
182  _eigenValues.clear();
183  _eigenValues.reserve(nImage);
184  for (int i = 0; i != nImage; ++i) {
185  _eigenValues.push_back(lambdaAndIndex[i].first);
186  }
187  //
188  // Contruct the first ncomp eigenimages in basis
189  //
190  int ncomp = 100; // number of components to keep
191 
192  _eigenImages.clear();
193  _eigenImages.reserve(ncomp < nImage ? ncomp : nImage);
194 
195  for (int i = 0; i < ncomp; ++i) {
196  if (i >= nImage) {
197  continue;
198  }
199 
200  int const ii = lambdaAndIndex[i].second; // the index after sorting (backwards) by eigenvalue
201 
202  std::shared_ptr<ImageT> eImage(new ImageT(_dimensions));
203  *eImage = static_cast<typename ImageT::Pixel>(0);
204 
205  for (int j = 0; j != nImage; ++j) {
206  int const jj = lambdaAndIndex[j].second; // the index after sorting (backwards) by eigenvalue
207  double const weight = Q(jj, ii) * (_constantWeight ? flux_bar / getFlux(jj) : 1);
208  eImage->scaledPlus(weight, *_imageList[jj]);
209  }
210  _eigenImages.push_back(eImage);
211  }
212 }
213 
214 namespace {
215 /*
216  * Fit a LinearCombinationKernel to an Image, allowing the coefficients of the components to vary
217  *
218  * return std::pair(best-fit kernel, std::pair(amp, chi^2))
219  */
220 template <typename MaskedImageT>
222  typename ImagePca<MaskedImageT>::ImageList const& eigenImages, // Eigen images
223  int nEigen, // Number of eigen images to use
224  MaskedImageT const& image // The image to be fit
225 ) {
226  typedef typename MaskedImageT::Image ImageT;
227 
228  if (nEigen == 0) {
229  throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "You must have at least one eigen image");
230  } else if (nEigen > static_cast<int>(eigenImages.size())) {
232  (boost::format("You only have %d eigen images (you asked for %d)") %
233  eigenImages.size() % nEigen)
234  .str());
235  }
236  /*
237  * Solve the linear problem image = sum x_i K_i + epsilon; we solve this for x_i by constructing the
238  * normal equations, A x = b
239  */
240  Eigen::MatrixXd A(nEigen, nEigen);
241  Eigen::VectorXd b(nEigen);
242 
243  for (int i = 0; i != nEigen; ++i) {
244  b(i) = innerProduct(*eigenImages[i]->getImage(), *image.getImage());
245 
246  for (int j = i; j != nEigen; ++j) {
247  A(i, j) = A(j, i) = innerProduct(*eigenImages[i]->getImage(), *eigenImages[j]->getImage());
248  }
249  }
250  Eigen::VectorXd x(nEigen);
251 
252  if (nEigen == 1) {
253  x(0) = b(0) / A(0, 0);
254  } else {
255  x = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
256  }
257  //
258  // Accumulate the best-fit-image in bestFitImage
259  //
260  std::shared_ptr<ImageT> bestFitImage = std::make_shared<ImageT>(eigenImages[0]->getDimensions());
261 
262  for (int i = 0; i != nEigen; ++i) {
263  bestFitImage->scaledPlus(x[i], *eigenImages[i]->getImage());
264  }
265 
266  return bestFitImage;
267 }
268 
269 template <typename ImageT>
270 double do_updateBadPixels(detail::basic_tag const&, typename ImagePca<ImageT>::ImageList const&,
271  std::vector<double> const&, typename ImagePca<ImageT>::ImageList const&,
272  unsigned long, int const) {
273  return 0.0;
274 }
275 
276 template <typename ImageT>
277 double do_updateBadPixels(detail::MaskedImage_tag const&,
278  typename ImagePca<ImageT>::ImageList const& imageList,
279  std::vector<double> const& fluxes, // fluxes of images
280  typename ImagePca<ImageT>::ImageList const& eigenImages, // Eigen images
281  unsigned long mask,
282  int const ncomp
283 ) {
284  int const nImage = imageList.size();
285  if (nImage == 0) {
287  "Please provide at least one Image for me to update");
288  }
289  lsst::geom::Extent2I dimensions = imageList[0]->getDimensions();
290  int const height = dimensions.getY();
291 
292  double maxChange = 0.0; // maximum change to the input images
293 
294  if (ncomp == 0) { // use mean of good pixels
295  typename ImageT::Image mean(dimensions); // desired mean image
296  image::Image<float> weight(mean.getDimensions()); // weight of each pixel
297 
298  for (int i = 0; i != nImage; ++i) {
299  double const flux_i = fluxes[i];
300 
301  for (int y = 0; y != height; ++y) {
302  typename ImageT::const_x_iterator iptr = imageList[i]->row_begin(y);
303  image::Image<float>::x_iterator wptr = weight.row_begin(y);
304  for (typename ImageT::Image::x_iterator mptr = mean.row_begin(y), end = mean.row_end(y);
305  mptr != end; ++mptr, ++iptr, ++wptr) {
306  if (!(iptr.mask() & mask) && iptr.variance() > 0.0) {
307  typename ImageT::Image::Pixel value = iptr.image() / flux_i;
308  float const var = iptr.variance() / (flux_i * flux_i);
309  float const ivar = 1.0 / var;
310 
311  if (std::isfinite(value * ivar)) {
312  *mptr += value * ivar;
313  *wptr += ivar;
314  }
315  }
316  }
317  }
318  }
319  //
320  // Calculate mean
321  //
322  for (int y = 0; y != height; ++y) {
323  image::Image<float>::x_iterator wptr = weight.row_begin(y);
324  for (typename ImageT::Image::x_iterator mptr = mean.row_begin(y), end = mean.row_end(y);
325  mptr != end; ++mptr, ++wptr) {
326  if (*wptr == 0.0) { // oops, no good values
327  ;
328  } else {
329  *mptr /= *wptr;
330  }
331  }
332  }
333  //
334  // Replace bad values by mean
335  //
336  for (int i = 0; i != nImage; ++i) {
337  double const flux_i = fluxes[i];
338 
339  for (int y = 0; y != height; ++y) {
340  typename ImageT::x_iterator iptr = imageList[i]->row_begin(y);
341  for (typename ImageT::Image::x_iterator mptr = mean.row_begin(y), end = mean.row_end(y);
342  mptr != end; ++mptr, ++iptr) {
343  if ((iptr.mask() & mask)) {
344  double const delta = ::fabs(flux_i * (*mptr) - iptr.image());
345  if (delta > maxChange) {
346  maxChange = delta;
347  }
348  iptr.image() = flux_i * (*mptr);
349  }
350  }
351  }
352  }
353  } else {
354  if (ncomp > static_cast<int>(eigenImages.size())) {
356  (boost::format("You only have %d eigen images (you asked for %d)") %
357  eigenImages.size() % ncomp)
358  .str());
359  }
360 
361  for (int i = 0; i != nImage; ++i) {
363  fitEigenImagesToImage(eigenImages, ncomp, *imageList[i]);
364 
365  for (int y = 0; y != height; ++y) {
366  typename ImageT::x_iterator iptr = imageList[i]->row_begin(y);
367  for (typename ImageT::Image::const_x_iterator fptr = fitted->row_begin(y),
368  end = fitted->row_end(y);
369  fptr != end; ++fptr, ++iptr) {
370  if (iptr.mask() & mask) {
371  double const delta = fabs(static_cast<double>(*fptr) - iptr.image());
372  if (delta > maxChange) {
373  maxChange = delta;
374  }
375 
376  iptr.image() = *fptr;
377  }
378  }
379  }
380  }
381  }
382 
383  return maxChange;
384 }
385 } // namespace
386 template <typename ImageT>
387 double ImagePca<ImageT>::updateBadPixels(unsigned long mask, int const ncomp) {
388  return do_updateBadPixels<ImageT>(typename ImageT::image_category(), _imageList, _fluxList, _eigenImages,
389  mask, ncomp);
390 }
391 
392 namespace {
393 template <typename T, typename U>
394 struct IsSame {
395  IsSame(T const&, U const&) {}
396  bool operator()() { return false; }
397 };
398 
399 template <typename T>
400 struct IsSame<T, T> {
401  IsSame(T const& im1, T const& im2) : _same(im1.row_begin(0) == im2.row_begin(0)) {}
402  bool operator()() { return _same; }
403 
404 private:
405  bool _same;
406 };
407 
408 // Test if two Images are identical; they need not be of the same type
409 template <typename Image1T, typename Image2T>
410 bool imagesAreIdentical(Image1T const& im1, Image2T const& im2) {
411  return IsSame<Image1T, Image2T>(im1, im2)();
412 }
413 } // namespace
414 template <typename Image1T, typename Image2T>
415 double innerProduct(Image1T const& lhs, Image2T const& rhs, int border) {
416  if (lhs.getWidth() <= 2 * border || lhs.getHeight() <= 2 * border) {
418  (boost::format("All image pixels are in the border of width %d: %dx%d") % border %
419  lhs.getWidth() % lhs.getHeight())
420  .str());
421  }
422 
423  double sum = 0.0;
424  //
425  // Handle I.I specially for efficiency, and to avoid advancing the iterator twice
426  //
427  if (imagesAreIdentical(lhs, rhs)) {
428  for (int y = border; y != lhs.getHeight() - border; ++y) {
429  for (typename Image1T::const_x_iterator lptr = lhs.row_begin(y) + border,
430  lend = lhs.row_end(y) - border;
431  lptr != lend; ++lptr) {
432  typename Image1T::Pixel val = *lptr;
433  if (std::isfinite(val)) {
434  sum += val * val;
435  }
436  }
437  }
438  } else {
439  if (lhs.getDimensions() != rhs.getDimensions()) {
441  (boost::format("Dimension mismatch: %dx%d v. %dx%d") % lhs.getWidth() %
442  lhs.getHeight() % rhs.getWidth() % rhs.getHeight())
443  .str());
444  }
445 
446  for (int y = border; y != lhs.getHeight() - border; ++y) {
447  typename Image2T::const_x_iterator rptr = rhs.row_begin(y) + border;
448  for (typename Image1T::const_x_iterator lptr = lhs.row_begin(y) + border,
449  lend = lhs.row_end(y) - border;
450  lptr != lend; ++lptr, ++rptr) {
451  double const tmp = (*lptr) * (*rptr);
452  if (std::isfinite(tmp)) {
453  sum += tmp;
454  }
455  }
456  }
457  }
458 
459  return sum;
460 }
461 
462 //
463 // Explicit instantiations
464 //
466 #define INSTANTIATE(T) \
467  template class ImagePca<Image<T> >; \
468  template double innerProduct(Image<T> const&, Image<T> const&, int); \
469  template class ImagePca<MaskedImage<T> >;
470 
471 #define INSTANTIATE2(T, U) \
472  template double innerProduct(Image<T> const&, Image<U> const&, int); \
473  template double innerProduct(Image<U> const&, Image<T> const&, int);
474 
477 INSTANTIATE(int)
478 INSTANTIATE(float)
479 INSTANTIATE(double)
480 
481 INSTANTIATE2(float, double) // the two types must be different
483 } // namespace image
484 } // namespace afw
485 } // namespace lsst
afw::table::PointKey< int > dimensions
Definition: GaussianPsf.cc:49
uint64_t * ptr
Definition: RangeSet.cc:88
T empty(T... args)
float Pixel
Typedefs to be used for pixel values.
Definition: common.h:37
ImageList getImageList() const
Return the list of images being analyzed.
Definition: ImagePca.cc:86
std::shared_ptr< ImageT > getMean() const
Return the mean of the images in ImagePca&#39;s list.
Definition: ImagePca.cc:91
double innerProduct(Image1T const &lhs, Image2T const &rhs, int const border=0)
Calculate the inner product of two images.
Definition: ImagePca.cc:415
table::Key< int > b
Reports attempts to exceed implementation-defined length limits for some classes. ...
Definition: Runtime.h:76
_view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: ImageBase.h:134
int y
Definition: SpanSet.cc:49
table::Key< int > a
STL namespace.
#define INSTANTIATE(TYPE)
Definition: ImagePca.cc:127
ImageT val
Definition: CR.cc:146
T end(T... args)
def dot(symb, c, r, frame=None, size=2, ctype=None, origin=afwImage.PARENT, args, kwargs)
Definition: ds9.py:101
ImagePca & operator=(ImagePca const &)
T push_back(T... args)
A base class for image defects.
A traits class for MaskedImage.
Definition: MaskedImage.h:51
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:168
T make_pair(T... args)
T isfinite(T... args)
double x
virtual void analyze()
Definition: ImagePca.cc:124
afw::table::Key< double > weight
ImagePca(bool constantWeight=true)
ctor
Definition: ImagePca.cc:47
T size(T... args)
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
Reports attempts to access elements outside a valid range of indices.
Definition: Runtime.h:89
afw::table::Key< afw::table::Array< MaskPixelT > > mask
lsst::geom::Extent2I const getDimensions() const
Return the dimension of the images being analyzed.
Definition: ImagePca.h:76
T begin(T... args)
void addImage(std::shared_ptr< ImageT > img, double flux=0.0)
Add an image to the set to be analyzed.
Definition: ImagePca.cc:65
T sort(T... args)
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
#define INSTANTIATE2(ImagePixelT1, ImagePixelT2)
Definition: MaskedImage.cc:689
int end
virtual double updateBadPixels(unsigned long mask, int const ncomp)
Update the bad pixels (i.e.
Definition: ImagePca.cc:387
T reserve(T... args)