LSST Applications g04c3c9f7ca+b6bf0e30ce,g2079a07aa2+3e9fd84d81,g20cdd03214+7686ae5cb8,g2305ad1205+9b0098d534,g2bbee38e9b+6c6beb4891,g337abbeb29+6c6beb4891,g3a166c0a6a+6c6beb4891,g4322eb9e3a+da50164803,g4e26c3aaff+ce16cedc45,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+9687a16dd4,g5c3423f6d4+7686ae5cb8,g607f77f49a+7686ae5cb8,g69c4c69b53+6c6beb4891,g858d7b2824+7686ae5cb8,g9963eaa53e+397fd696c8,g99cad8db69+2f9f763ee9,g9ddcbc5298+9a081db1e4,ga1e77700b3+2cbb763275,ga8d6e2f9c9+cbb649eb12,gadfd92a7e4+aec2f3b930,gae0086650b+585e252eca,gb0e22166c9+0e73c8378f,gbb8dafda3b+024d427455,gc120e1dc64+fcd4152152,gc28159a63d+6c6beb4891,gcdd4ae20e8+3b72bea61f,gcf0d15dbbd+3b72bea61f,gd1535ee943+38e04e5621,gd71bba6883+4b9ebe24e7,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+6c6beb4891,gfbcc870c63+b310236976,w.2024.24
LSST Data Management Base Package
Loading...
Searching...
No Matches
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/Eigenvalues"
35
38
39namespace afwMath = lsst::afw::math;
40
41namespace lsst {
42namespace afw {
43namespace image {
44
45template <typename ImageT>
46ImagePca<ImageT>::ImagePca(bool constantWeight)
47 : _imageList(),
48 _fluxList(),
49 _dimensions(0, 0),
50 _constantWeight(constantWeight),
51 _eigenValues(std::vector<double>()),
52 _eigenImages(ImageList()) {}
53
54template <typename ImageT>
55ImagePca<ImageT>::ImagePca(ImagePca const&) = default;
56template <typename ImageT>
58template <typename ImageT>
60template <typename ImageT>
62
63template <typename ImageT>
65 if (_imageList.empty()) {
66 _dimensions = img->getDimensions();
67 } else {
68 if (getDimensions() != img->getDimensions()) {
70 (boost::format("Dimension mismatch: saw %dx%d; expected %dx%d") %
71 img->getWidth() % img->getHeight() % _dimensions.getX() % _dimensions.getY())
72 .str());
73 }
74 }
75
76 if (flux == 0.0) {
77 throw LSST_EXCEPT(lsst::pex::exceptions::OutOfRangeError, "Flux may not be zero");
78 }
79
80 _imageList.push_back(img);
81 _fluxList.push_back(flux);
82}
83
84template <typename ImageT>
86 return _imageList;
87}
88
89template <typename ImageT>
91 if (_imageList.empty()) {
92 throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "You haven't provided any images");
93 }
94
95 std::shared_ptr<ImageT> mean(new ImageT(getDimensions()));
96 *mean = static_cast<typename ImageT::Pixel>(0);
97
98 for (typename ImageList::const_iterator ptr = _imageList.begin(), end = _imageList.end(); ptr != end;
99 ++ptr) {
100 *mean += **ptr;
101 }
102 *mean /= _imageList.size();
103
104 return mean;
105}
106
107namespace {
108/*
109 * Analyze the images in an ImagePca, calculating the PCA decomposition (== Karhunen-Lo\`eve basis)
110 *
111 * The notation is that in chapter 7 of Gyula Szokoly's thesis at JHU
112 */
113template <typename T>
114struct SortEvalueDecreasing {
115 bool operator()(std::pair<T, int> const& a, std::pair<T, int> const& b) {
116 return a.first > b.first; // N.b. sort on greater
117 }
118};
119} // namespace
120
121template <typename ImageT>
123 int const nImage = _imageList.size();
124
125 if (nImage == 0) {
126 throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "No images provided for PCA analysis");
127 }
128 /*
129 * Eigen doesn't like 1x1 matrices, but we don't really need it to handle a single matrix...
130 */
131 if (nImage == 1) {
132 _eigenImages.clear();
133 _eigenImages.push_back(std::shared_ptr<ImageT>(new ImageT(*_imageList[0], true)));
134
135 _eigenValues.clear();
136 _eigenValues.push_back(1.0);
137
138 return;
139 }
140 /*
141 * Find the eigenvectors/values of the scalar product matrix, R' (Eq. 7.4)
142 */
143 Eigen::MatrixXd R(nImage, nImage); // residuals' inner products
144
145 double flux_bar = 0; // mean of flux for all regions
146 for (int i = 0; i != nImage; ++i) {
147 typename GetImage<ImageT>::type const& im_i = *GetImage<ImageT>::getImage(_imageList[i]);
148 double const flux_i = getFlux(i);
149 flux_bar += flux_i;
150
151 for (int j = i; j != nImage; ++j) {
152 typename GetImage<ImageT>::type const& im_j = *GetImage<ImageT>::getImage(_imageList[j]);
153 double const flux_j = getFlux(j);
154
155 double dot = innerProduct(im_i, im_j);
156 if (_constantWeight) {
157 dot /= flux_i * flux_j;
158 }
159 R(i, j) = R(j, i) = dot / nImage;
160 }
161 }
162 flux_bar /= nImage;
163 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(R);
164 Eigen::MatrixXd const& Q = eVecValues.eigenvectors();
165 Eigen::VectorXd const& lambda = eVecValues.eigenvalues();
166 //
167 // We need to sort the eigenValues, and remember the permutation we applied to the eigenImages
168 // We'll use the vector lambdaAndIndex to achieve this
169 //
170 std::vector<std::pair<double, int> > lambdaAndIndex; // pairs (eValue, index)
171 lambdaAndIndex.reserve(nImage);
172
173 for (int i = 0; i != nImage; ++i) {
174 lambdaAndIndex.emplace_back(lambda(i), i);
175 }
176 std::sort(lambdaAndIndex.begin(), lambdaAndIndex.end(), SortEvalueDecreasing<double>());
177 //
178 // Save the (sorted) eigen values
179 //
180 _eigenValues.clear();
181 _eigenValues.reserve(nImage);
182 for (int i = 0; i != nImage; ++i) {
183 _eigenValues.push_back(lambdaAndIndex[i].first);
184 }
185 //
186 // Contruct the first ncomp eigenimages in basis
187 //
188 int ncomp = 100; // number of components to keep
189
190 _eigenImages.clear();
191 _eigenImages.reserve(ncomp < nImage ? ncomp : nImage);
192
193 for (int i = 0; i < ncomp; ++i) {
194 if (i >= nImage) {
195 continue;
196 }
197
198 int const ii = lambdaAndIndex[i].second; // the index after sorting (backwards) by eigenvalue
199
200 std::shared_ptr<ImageT> eImage(new ImageT(_dimensions));
201 *eImage = static_cast<typename ImageT::Pixel>(0);
202
203 for (int j = 0; j != nImage; ++j) {
204 int const jj = lambdaAndIndex[j].second; // the index after sorting (backwards) by eigenvalue
205 double const weight = Q(jj, ii) * (_constantWeight ? flux_bar / getFlux(jj) : 1);
206 eImage->scaledPlus(weight, *_imageList[jj]);
207 }
208 _eigenImages.push_back(eImage);
209 }
210}
211
212namespace {
213/*
214 * Fit a LinearCombinationKernel to an Image, allowing the coefficients of the components to vary
215 *
216 * return std::pair(best-fit kernel, std::pair(amp, chi^2))
217 */
218template <typename MaskedImageT>
220 typename ImagePca<MaskedImageT>::ImageList const& eigenImages, // Eigen images
221 int nEigen, // Number of eigen images to use
222 MaskedImageT const& image // The image to be fit
223) {
224 using ImageT = typename MaskedImageT::Image;
225
226 if (nEigen == 0) {
227 throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "You must have at least one eigen image");
228 } else if (nEigen > static_cast<int>(eigenImages.size())) {
230 (boost::format("You only have %d eigen images (you asked for %d)") %
231 eigenImages.size() % nEigen)
232 .str());
233 }
234 /*
235 * Solve the linear problem image = sum x_i K_i + epsilon; we solve this for x_i by constructing the
236 * normal equations, A x = b
237 */
238 Eigen::MatrixXd A(nEigen, nEigen);
239 Eigen::VectorXd b(nEigen);
240
241 for (int i = 0; i != nEigen; ++i) {
242 b(i) = innerProduct(*eigenImages[i]->getImage(), *image.getImage());
243
244 for (int j = i; j != nEigen; ++j) {
245 A(i, j) = A(j, i) = innerProduct(*eigenImages[i]->getImage(), *eigenImages[j]->getImage());
246 }
247 }
248 Eigen::VectorXd x(nEigen);
249
250 if (nEigen == 1) {
251 x(0) = b(0) / A(0, 0);
252 } else {
253 x = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
254 }
255 //
256 // Accumulate the best-fit-image in bestFitImage
257 //
258 std::shared_ptr<ImageT> bestFitImage = std::make_shared<ImageT>(eigenImages[0]->getDimensions());
259
260 for (int i = 0; i != nEigen; ++i) {
261 bestFitImage->scaledPlus(x[i], *eigenImages[i]->getImage());
262 }
263
264 return bestFitImage;
265}
266
267template <typename ImageT>
268double do_updateBadPixels(detail::basic_tag const&, typename ImagePca<ImageT>::ImageList const&,
269 std::vector<double> const&, typename ImagePca<ImageT>::ImageList const&,
270 unsigned long, int const) {
271 return 0.0;
272}
273
274template <typename ImageT>
275double do_updateBadPixels(detail::MaskedImage_tag const&,
276 typename ImagePca<ImageT>::ImageList const& imageList,
277 std::vector<double> const& fluxes, // fluxes of images
278 typename ImagePca<ImageT>::ImageList const& eigenImages, // Eigen images
279 unsigned long mask,
280 int const ncomp
281) {
282 int const nImage = imageList.size();
283 if (nImage == 0) {
285 "Please provide at least one Image for me to update");
286 }
287 lsst::geom::Extent2I dimensions = imageList[0]->getDimensions();
288 int const height = dimensions.getY();
289
290 double maxChange = 0.0; // maximum change to the input images
291
292 if (ncomp == 0) { // use mean of good pixels
293 typename ImageT::Image mean(dimensions); // desired mean image
294 image::Image<float> weight(mean.getDimensions()); // weight of each pixel
295
296 for (int i = 0; i != nImage; ++i) {
297 double const flux_i = fluxes[i];
298
299 for (int y = 0; y != height; ++y) {
300 typename ImageT::const_x_iterator iptr = imageList[i]->row_begin(y);
302 for (typename ImageT::Image::x_iterator mptr = mean.row_begin(y), end = mean.row_end(y);
303 mptr != end; ++mptr, ++iptr, ++wptr) {
304 if (!(iptr.mask() & mask) && iptr.variance() > 0.0) {
305 typename ImageT::Image::Pixel value = iptr.image() / flux_i;
306 float const var = iptr.variance() / (flux_i * flux_i);
307 float const ivar = 1.0 / var;
308
309 if (std::isfinite(value * ivar)) {
310 *mptr += value * ivar;
311 *wptr += ivar;
312 }
313 }
314 }
315 }
316 }
317 //
318 // Calculate mean
319 //
320 for (int y = 0; y != height; ++y) {
322 for (typename ImageT::Image::x_iterator mptr = mean.row_begin(y), end = mean.row_end(y);
323 mptr != end; ++mptr, ++wptr) {
324 if (*wptr == 0.0) { // oops, no good values
325 ;
326 } else {
327 *mptr /= *wptr;
328 }
329 }
330 }
331 //
332 // Replace bad values by mean
333 //
334 for (int i = 0; i != nImage; ++i) {
335 double const flux_i = fluxes[i];
336
337 for (int y = 0; y != height; ++y) {
338 typename ImageT::x_iterator iptr = imageList[i]->row_begin(y);
339 for (typename ImageT::Image::x_iterator mptr = mean.row_begin(y), end = mean.row_end(y);
340 mptr != end; ++mptr, ++iptr) {
341 if ((iptr.mask() & mask)) {
342 double const delta = ::fabs(flux_i * (*mptr) - iptr.image());
343 if (delta > maxChange) {
344 maxChange = delta;
345 }
346 iptr.image() = flux_i * (*mptr);
347 }
348 }
349 }
350 }
351 } else {
352 if (ncomp > static_cast<int>(eigenImages.size())) {
354 (boost::format("You only have %d eigen images (you asked for %d)") %
355 eigenImages.size() % ncomp)
356 .str());
357 }
358
359 for (int i = 0; i != nImage; ++i) {
361 fitEigenImagesToImage(eigenImages, ncomp, *imageList[i]);
362
363 for (int y = 0; y != height; ++y) {
364 typename ImageT::x_iterator iptr = imageList[i]->row_begin(y);
365 for (typename ImageT::Image::const_x_iterator fptr = fitted->row_begin(y),
366 end = fitted->row_end(y);
367 fptr != end; ++fptr, ++iptr) {
368 if (iptr.mask() & mask) {
369 double const delta = fabs(static_cast<double>(*fptr) - iptr.image());
370 if (delta > maxChange) {
371 maxChange = delta;
372 }
373
374 iptr.image() = *fptr;
375 }
376 }
377 }
378 }
379 }
380
381 return maxChange;
382}
383} // namespace
384template <typename ImageT>
385double ImagePca<ImageT>::updateBadPixels(unsigned long mask, int const ncomp) {
386 return do_updateBadPixels<ImageT>(typename ImageT::image_category(), _imageList, _fluxList, _eigenImages,
387 mask, ncomp);
388}
389
390namespace {
391template <typename T, typename U>
392struct IsSame {
393 IsSame(T const&, U const&) {}
394 bool operator()() { return false; }
395};
396
397template <typename T>
398struct IsSame<T, T> {
399 IsSame(T const& im1, T const& im2) : _same(im1.row_begin(0) == im2.row_begin(0)) {}
400 bool operator()() { return _same; }
401
402private:
403 bool _same;
404};
405
406// Test if two Images are identical; they need not be of the same type
407template <typename Image1T, typename Image2T>
408bool imagesAreIdentical(Image1T const& im1, Image2T const& im2) {
409 return IsSame<Image1T, Image2T>(im1, im2)();
410}
411} // namespace
412template <typename Image1T, typename Image2T>
413double innerProduct(Image1T const& lhs, Image2T const& rhs, int border) {
414 if (lhs.getWidth() <= 2 * border || lhs.getHeight() <= 2 * border) {
416 (boost::format("All image pixels are in the border of width %d: %dx%d") % border %
417 lhs.getWidth() % lhs.getHeight())
418 .str());
419 }
420
421 double sum = 0.0;
422 //
423 // Handle I.I specially for efficiency, and to avoid advancing the iterator twice
424 //
425 if (imagesAreIdentical(lhs, rhs)) {
426 for (int y = border; y != lhs.getHeight() - border; ++y) {
427 for (typename Image1T::const_x_iterator lptr = lhs.row_begin(y) + border,
428 lend = lhs.row_end(y) - border;
429 lptr != lend; ++lptr) {
430 typename Image1T::Pixel val = *lptr;
431 if (std::isfinite(val)) {
432 sum += val * val;
433 }
434 }
435 }
436 } else {
437 if (lhs.getDimensions() != rhs.getDimensions()) {
439 (boost::format("Dimension mismatch: %dx%d v. %dx%d") % lhs.getWidth() %
440 lhs.getHeight() % rhs.getWidth() % rhs.getHeight())
441 .str());
442 }
443
444 for (int y = border; y != lhs.getHeight() - border; ++y) {
445 typename Image2T::const_x_iterator rptr = rhs.row_begin(y) + border;
446 for (typename Image1T::const_x_iterator lptr = lhs.row_begin(y) + border,
447 lend = lhs.row_end(y) - border;
448 lptr != lend; ++lptr, ++rptr) {
449 double const tmp = (*lptr) * (*rptr);
450 if (std::isfinite(tmp)) {
451 sum += tmp;
452 }
453 }
454 }
455 }
456
457 return sum;
458}
459
460//
461// Explicit instantiations
462//
464#define INSTANTIATE(T) \
465 template class ImagePca<Image<T> >; \
466 template double innerProduct(Image<T> const&, Image<T> const&, int); \
467 template class ImagePca<MaskedImage<T> >;
468
469#define INSTANTIATE2(T, U) \
470 template double innerProduct(Image<T> const&, Image<U> const&, int); \
471 template double innerProduct(Image<U> const&, Image<T> const&, int);
472
475INSTANTIATE(int)
476INSTANTIATE(float)
477INSTANTIATE(double)
478
479INSTANTIATE2(float, double) // the two types must be different
481} // namespace image
482} // namespace afw
483} // namespace lsst
int end
#define INSTANTIATE(FROMSYS, TOSYS)
Definition Detector.cc:509
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
afw::table::Key< afw::table::Array< MaskPixelT > > mask
#define INSTANTIATE2(ImagePixelT1, ImagePixelT2)
uint64_t * ptr
Definition RangeSet.cc:95
int y
Definition SpanSet.cc:48
table::Key< int > b
T begin(T... args)
x_iterator row_begin(int y) const
Return an x_iterator to the start of the y'th row.
Definition ImageBase.h:401
A class to represent a 2-dimensional array of pixels.
Definition Image.h:51
virtual double updateBadPixels(unsigned long mask, int const ncomp)
Update the bad pixels (i.e.
Definition ImagePca.cc:385
void addImage(std::shared_ptr< ImageT > img, double flux=0.0)
Add an image to the set to be analyzed.
Definition ImagePca.cc:64
ImageList getImageList() const
Return the list of images being analyzed.
Definition ImagePca.cc:85
ImagePca(bool constantWeight=true)
ctor
Definition ImagePca.cc:46
std::shared_ptr< ImageT > getMean() const
Return the mean of the images in ImagePca's list.
Definition ImagePca.cc:90
std::vector< std::shared_ptr< ImageT > > ImageList
Definition ImagePca.h:47
ImagePca & operator=(ImagePca const &)
Reports attempts to exceed implementation-defined length limits for some classes.
Definition Runtime.h:76
Reports attempts to access elements outside a valid range of indices.
Definition Runtime.h:89
T emplace_back(T... args)
T end(T... args)
T fabs(T... args)
T isfinite(T... args)
double innerProduct(Image1T const &lhs, Image2T const &rhs, int const border=0)
Calculate the inner product of two images.
Definition ImagePca.cc:413
STL namespace.
T reserve(T... args)
T size(T... args)
T sort(T... args)
afw::table::Key< double > weight
ImageT val
Definition CR.cc:146