LSST Applications g0265f82a02+c6dfa2ddaf,g1162b98a3f+b2075782a9,g2079a07aa2+1b2e822518,g2bbee38e9b+c6dfa2ddaf,g337abbeb29+c6dfa2ddaf,g3ddfee87b4+a60788ef87,g50ff169b8f+2eb0e556e8,g52b1c1532d+90ebb246c7,g555ede804d+a60788ef87,g591dd9f2cf+ba8caea58f,g5ec818987f+864ee9cddb,g858d7b2824+9ee1ab4172,g876c692160+a40945ebb7,g8a8a8dda67+90ebb246c7,g8cdfe0ae6a+4fd9e222a8,g99cad8db69+5e309b7bc6,g9ddcbc5298+a1346535a5,ga1e77700b3+df8f93165b,ga8c6da7877+aa12a14d27,gae46bcf261+c6dfa2ddaf,gb0e22166c9+8634eb87fb,gb3f2274832+d0da15e3be,gba4ed39666+1ac82b564f,gbb8dafda3b+5dfd9c994b,gbeb006f7da+97157f9740,gc28159a63d+c6dfa2ddaf,gc86a011abf+9ee1ab4172,gcf0d15dbbd+a60788ef87,gdaeeff99f8+1cafcb7cd4,gdc0c513512+9ee1ab4172,ge79ae78c31+c6dfa2ddaf,geb67518f79+ba1859f325,geb961e4c1e+f9439d1e6f,gee10cc3b42+90ebb246c7,gf1cff7945b+9ee1ab4172,w.2024.12
LSST Data Management Base Package
Loading...
Searching...
No Matches
Public Types | Public Member Functions | List of all members
lsst::afw::image::ImagePca< ImageT > Class Template Reference

#include <ImagePca.h>

Inheritance diagram for lsst::afw::image::ImagePca< ImageT >:
lsst::ip::diffim::detail::KernelPca< Image< lsst::afw::math::Kernel::Pixel > > lsst::ip::diffim::detail::KernelPca< ImageT > lsst::meas::algorithms::PsfImagePca< ImageT >

Public Types

using ImageList = std::vector<std::shared_ptr<ImageT>>
 

Public Member Functions

 ImagePca (bool constantWeight=true)
 ctor
 
virtual ~ImagePca ()=default
 
 ImagePca (ImagePca const &)
 
 ImagePca (ImagePca &&)
 
ImagePcaoperator= (ImagePca const &)
 
ImagePcaoperator= (ImagePca &&)
 
void addImage (std::shared_ptr< ImageT > img, double flux=0.0)
 Add an image to the set to be analyzed.
 
ImageList getImageList () const
 Return the list of images being analyzed.
 
lsst::geom::Extent2I const getDimensions () const
 Return the dimension of the images being analyzed.
 
std::shared_ptr< ImageT > getMean () const
 Return the mean of the images in ImagePca's list.
 
virtual void analyze ()
 
virtual double updateBadPixels (unsigned long mask, int const ncomp)
 Update the bad pixels (i.e.
 
std::vector< double > const & getEigenValues () const
 Return Eigen values.
 
ImageList const & getEigenImages () const
 Return Eigen images.
 

Detailed Description

template<typename ImageT>
class lsst::afw::image::ImagePca< ImageT >

Definition at line 45 of file ImagePca.h.

Member Typedef Documentation

◆ ImageList

template<typename ImageT >
using lsst::afw::image::ImagePca< ImageT >::ImageList = std::vector<std::shared_ptr<ImageT>>

Definition at line 47 of file ImagePca.h.

Constructor & Destructor Documentation

◆ ImagePca() [1/3]

template<typename ImageT >
lsst::afw::image::ImagePca< ImageT >::ImagePca ( bool constantWeight = true)
explicit

ctor

Parameters
constantWeightShould all stars be weighted equally?

Definition at line 46 of file ImagePca.cc.

47 : _imageList(),
48 _fluxList(),
49 _dimensions(0, 0),
50 _constantWeight(constantWeight),
51 _eigenValues(std::vector<double>()),
52 _eigenImages(ImageList()) {}
std::vector< std::shared_ptr< ImageT > > ImageList
Definition ImagePca.h:47

◆ ~ImagePca()

template<typename ImageT >
virtual lsst::afw::image::ImagePca< ImageT >::~ImagePca ( )
virtualdefault

◆ ImagePca() [2/3]

template<typename ImageT >
lsst::afw::image::ImagePca< ImageT >::ImagePca ( ImagePca< ImageT > const & )
default

◆ ImagePca() [3/3]

template<typename ImageT >
lsst::afw::image::ImagePca< ImageT >::ImagePca ( ImagePca< ImageT > && )
default

Member Function Documentation

◆ addImage()

template<typename ImageT >
void lsst::afw::image::ImagePca< ImageT >::addImage ( std::shared_ptr< ImageT > img,
double flux = 0.0 )

Add an image to the set to be analyzed.

Parameters
imgImage to add to set
fluxImage's flux
Exceptions
lsst::pex::exceptions::LengthErrorif all the images aren't the same size

Definition at line 64 of file ImagePca.cc.

64 {
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}
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
lsst::geom::Extent2I const getDimensions() const
Return the dimension of the images being analyzed.
Definition ImagePca.h:75
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 empty(T... args)
T push_back(T... args)

◆ analyze()

template<typename ImageT >
void lsst::afw::image::ImagePca< ImageT >::analyze ( )
virtual

Reimplemented in lsst::ip::diffim::detail::KernelPca< Image< lsst::afw::math::Kernel::Pixel > >, lsst::ip::diffim::detail::KernelPca< ImageT >, and lsst::meas::algorithms::PsfImagePca< ImageT >.

Definition at line 122 of file ImagePca.cc.

122 {
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}
T begin(T... args)
T clear(T... args)
T emplace_back(T... args)
T end(T... args)
dot(symb, c, r, frame=None, size=2, ctype=None, origin=afwImage.PARENT, *args, **kwargs)
Definition ds9.py:100
double innerProduct(Image1T const &lhs, Image2T const &rhs, int const border=0)
Calculate the inner product of two images.
Definition ImagePca.cc:413
T reserve(T... args)
T size(T... args)
T sort(T... args)
afw::table::Key< double > weight

◆ getDimensions()

template<typename ImageT >
lsst::geom::Extent2I const lsst::afw::image::ImagePca< ImageT >::getDimensions ( ) const
inline

Return the dimension of the images being analyzed.

Definition at line 75 of file ImagePca.h.

75{ return _dimensions; }

◆ getEigenImages()

template<typename ImageT >
ImageList const & lsst::afw::image::ImagePca< ImageT >::getEigenImages ( ) const
inline

Return Eigen images.

Definition at line 100 of file ImagePca.h.

100{ return _eigenImages; }

◆ getEigenValues()

template<typename ImageT >
std::vector< double > const & lsst::afw::image::ImagePca< ImageT >::getEigenValues ( ) const
inline

Return Eigen values.

Definition at line 98 of file ImagePca.h.

98{ return _eigenValues; }

◆ getImageList()

template<typename ImageT >
ImagePca< ImageT >::ImageList lsst::afw::image::ImagePca< ImageT >::getImageList ( ) const

Return the list of images being analyzed.

Definition at line 85 of file ImagePca.cc.

85 {
86 return _imageList;
87}

◆ getMean()

template<typename ImageT >
std::shared_ptr< ImageT > lsst::afw::image::ImagePca< ImageT >::getMean ( ) const

Return the mean of the images in ImagePca's list.

Definition at line 90 of file ImagePca.cc.

90 {
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}
int end
uint64_t * ptr
Definition RangeSet.cc:95

◆ operator=() [1/2]

template<typename ImageT >
ImagePca< ImageT > & lsst::afw::image::ImagePca< ImageT >::operator= ( ImagePca< ImageT > && )
default

◆ operator=() [2/2]

template<typename ImageT >
ImagePca< ImageT > & lsst::afw::image::ImagePca< ImageT >::operator= ( ImagePca< ImageT > const & )
default

◆ updateBadPixels()

template<typename ImageT >
double lsst::afw::image::ImagePca< ImageT >::updateBadPixels ( unsigned long mask,
int const ncomp )
virtual

Update the bad pixels (i.e.

those for which (value & mask) != 0) based on the current PCA decomposition; if none is available, use the mean of the good pixels

Parameters
maskMask defining bad pixels
ncompNumber of components to use in estimate
Returns
the maximum change made to any pixel

N.b. the work is actually done in do_updateBadPixels as the code only makes sense and compiles when we are doing a PCA on a set of MaskedImages

Definition at line 385 of file ImagePca.cc.

385 {
386 return do_updateBadPixels<ImageT>(typename ImageT::image_category(), _imageList, _fluxList, _eigenImages,
387 mask, ncomp);
388}
afw::table::Key< afw::table::Array< MaskPixelT > > mask

The documentation for this class was generated from the following files: