LSSTApplications  16.0-1-gce273f5+20,16.0-10-g1758552,16.0-10-gc1446dd+21,16.0-12-g569485f+3,16.0-13-g5f20d24+1,16.0-13-g80874fd+2,16.0-13-gb122224+12,16.0-14-g08f9460+4,16.0-14-g8a3b804,16.0-14-ga5060d2,16.0-15-g77ef378+7,16.0-18-gdf247dd+2,16.0-18-ge18fa5b,16.0-18-ge4151178,16.0-2-g0febb12+16,16.0-2-g9d5294e+46,16.0-2-gc6e0ed0+5,16.0-23-ge8a9b866+3,16.0-3-g404ea43+13,16.0-3-gbc759ec+19,16.0-3-gcfd6c53+44,16.0-3-ge00e371,16.0-4-g03cf288+35,16.0-4-g13a27c5+21,16.0-4-g5f3a788+15,16.0-4-ga3eb747+5,16.0-5-g1991253+21,16.0-5-g1e9226d+3,16.0-5-g6a53317,16.0-5-g865efd9+23,16.0-5-gb3f8a4b+53,16.0-5-gd0f1235+10,16.0-52-gad2f36c79,16.0-7-g6043bfc+9,16.0-7-gd2eeba5+3,16.0-7-gde5bd64+3,16.0-8-g0e813a6+1,16.0-9-g52c50f7,16.0-9-g73415e6,master-g5768c874b9+5,w.2018.41
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 
33 #include "lsst/afw.h"
35 
36 namespace lsst {
37 namespace meas {
38 namespace algorithms {
39 
40 template <typename ImageT>
42  Super::analyze();
43 
44  typename Super::ImageList const &eImageList = this->getEigenImages();
45  typename Super::ImageList::const_iterator iter = eImageList.begin(), end = eImageList.end();
46  for (size_t i = 0; iter != end; ++i, ++iter) {
47  PTR(ImageT) eImage = *iter;
48 
49  /*
50  * Normalise eigenImages to have a maximum of 1.0. For n > 0 they
51  * (should) have mean == 0, so we can't use that to normalize
52  */
54  double const min = stats.getValue(afw::math::MIN);
55  double const max = stats.getValue(afw::math::MAX);
56 
57  double const extreme = (fabs(min) > max) ? min : max;
58  if (extreme != 0.0) {
59  *eImage /= extreme;
60  }
61 
62  /*
63  * Estimate and subtract the mean background level from the i > 0
64  * eigen images; if we don't do that then PSF variation can get mixed
65  * with subtle variations in the background and potentially amplify
66  * them disasterously.
67  *
68  * It is not at all clear that doing this is a good idea; it'd be
69  * better to get the sky level right in the first place.
70  */
71  if (i > 0 && _border > 0) { /* not the zeroth KL component */
72  int const height = eImage->getHeight();
73  int const width = eImage->getWidth();
74  double background;
75  if (2 * _border >= std::min(height, width)) {
76  // _Border consumes the entire image
79  .getValue();
80  } else {
81  // Use the median of the edge pixels
82 
83  // If ImageT is a MaskedImage, unpack the Image
86 
87  int const nEdge = width * height - (width - 2 * _border) * (height - 2 * _border);
88  std::vector<double> edgePixels(nEdge);
89 
90  std::vector<double>::iterator bi = edgePixels.begin();
91 
92  typedef typename afw::image::GetImage<ImageT>::type::x_iterator imIter;
93  int y = 0;
94  for (; y != _border; ++y) { // Bottom border of eImage
95  for (imIter ptr = eImageIm->row_begin(y), end = eImageIm->row_end(y); ptr != end;
96  ++ptr, ++bi) {
97  *bi = *ptr;
98  }
99  }
100  for (; y != height - _border; ++y) { // Left and right borders of eImage
101  for (imIter ptr = eImageIm->row_begin(y), end = eImageIm->x_at(_border, y); ptr != end;
102  ++ptr, ++bi) {
103  *bi = *ptr;
104  }
105  for (imIter ptr = eImageIm->x_at(width - _border, y), end = eImageIm->row_end(y);
106  ptr != end; ++ptr, ++bi) {
107  *bi = *ptr;
108  }
109  }
110  for (; y != height; ++y) { // Top border of eImage
111  for (imIter ptr = eImageIm->row_begin(y), end = eImageIm->row_end(y); ptr != end;
112  ++ptr, ++bi) {
113  *bi = *ptr;
114  }
115  }
116  assert(std::distance(edgePixels.begin(), bi) == nEdge);
117 
118  background = afw::math::makeStatistics(edgePixels, afw::math::MEDIAN).getValue();
119  }
120  *eImage -= background;
121  }
122  }
123 }
124 
125 #define INSTANTIATE_IMAGE(IMAGE) template class PsfImagePca<IMAGE>;
126 
127 #define INSTANTIATE(TYPE) \
128  INSTANTIATE_IMAGE(afw::image::Image<TYPE>); \
129  INSTANTIATE_IMAGE(afw::image::MaskedImage<TYPE>);
130 
131 INSTANTIATE(float);
132 
133 } // namespace algorithms
134 } // namespace meas
135 } // namespace lsst
uint64_t * ptr
Definition: RangeSet.cc:88
T distance(T... args)
estimate sample minimum
Definition: Statistics.h:76
int y
Definition: SpanSet.cc:49
virtual void analyze()
Generate eigenimages that are normalised and background-subtracted.
Definition: ImagePca.cc:41
estimate sample maximum
Definition: Statistics.h:77
T end(T... args)
int min
A class to evaluate image statistics.
Definition: Statistics.h:215
#define PTR(...)
Definition: base.h:41
T min(T... args)
estimate sample median
Definition: Statistics.h:70
A base class for image defects.
Definition: cameraGeom.dox:3
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
#define INSTANTIATE(TYPE)
Definition: ImagePca.cc:127
Class for doing PCA on PSF stars.
int max
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1056
T begin(T... args)
int end