lsst::afw::image::ImagePca< ImageT > Class Template Reference

#include <ImagePca.h>

Inheritance diagram for lsst::afw::image::ImagePca< ImageT >:

Public Types
typedef std::vector< std::shared_ptr< ImageT > >	ImageList

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

Detailed Description

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

Definition at line 46 of file ImagePca.h.

Member Typedef Documentation

ImageList

template<typename ImageT >

typedef std::vector<std::shared_ptr<ImageT> > lsst::afw::image::ImagePca< ImageT >::ImageList

Definition at line 48 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

constantWeight

Should all stars be weighted equally?

Definition at line 47 of file ImagePca.cc.

        : _imageList(),
          _fluxList(),
          _dimensions(0, 0),
          _constantWeight(constantWeight),
          _eigenValues(std::vector<double>()),
          _eigenImages(ImageList()) {}

~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

img	Image to add to set
flux	Image's flux

Exceptions

lsst::pex::exceptions::LengthError

if all the images aren't the same size

Definition at line 65 of file ImagePca.cc.

                                                                      {
    if (_imageList.empty()) {
        _dimensions = img->getDimensions();
    } else {
        if (getDimensions() != img->getDimensions()) {
            throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                              (boost::format("Dimension mismatch: saw %dx%d; expected %dx%d") %
                               img->getWidth() % img->getHeight() % _dimensions.getX() % _dimensions.getY())
                                      .str());
        }
    }
 
    if (flux == 0.0) {
        throw LSST_EXCEPT(lsst::pex::exceptions::OutOfRangeError, "Flux may not be zero");
    }
 
    _imageList.push_back(img);
    _fluxList.push_back(flux);
}

analyze()

template<typename ImageT >

void lsst::afw::image::ImagePca< ImageT >::analyze

virtual

Reimplemented in lsst::meas::algorithms::PsfImagePca< ImageT >, and lsst::ip::diffim::detail::KernelPca< ImageT >.

Definition at line 124 of file ImagePca.cc.

                               {
    int const nImage = _imageList.size();
 
    if (nImage == 0) {
        throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "No images provided for PCA analysis");
    }
    /*
     * Eigen doesn't like 1x1 matrices, but we don't really need it to handle a single matrix...
     */
    if (nImage == 1) {
        _eigenImages.clear();
        _eigenImages.push_back(std::shared_ptr<ImageT>(new ImageT(*_imageList[0], true)));
 
        _eigenValues.clear();
        _eigenValues.push_back(1.0);
 
        return;
    }
    /*
     * Find the eigenvectors/values of the scalar product matrix, R' (Eq. 7.4)
     */
    Eigen::MatrixXd R(nImage, nImage);  // residuals' inner products
 
    double flux_bar = 0;  // mean of flux for all regions
    for (int i = 0; i != nImage; ++i) {
        typename GetImage<ImageT>::type const& im_i = *GetImage<ImageT>::getImage(_imageList[i]);
        double const flux_i = getFlux(i);
        flux_bar += flux_i;
 
        for (int j = i; j != nImage; ++j) {
            typename GetImage<ImageT>::type const& im_j = *GetImage<ImageT>::getImage(_imageList[j]);
            double const flux_j = getFlux(j);
 
            double dot = innerProduct(im_i, im_j);
            if (_constantWeight) {
                dot /= flux_i * flux_j;
            }
            R(i, j) = R(j, i) = dot / nImage;
        }
    }
    flux_bar /= nImage;
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eVecValues(R);
    Eigen::MatrixXd const& Q = eVecValues.eigenvectors();
    Eigen::VectorXd const& lambda = eVecValues.eigenvalues();
    //
    // We need to sort the eigenValues, and remember the permutation we applied to the eigenImages
    // We'll use the vector lambdaAndIndex to achieve this
    //
    std::vector<std::pair<double, int> > lambdaAndIndex;  // pairs (eValue, index)
    lambdaAndIndex.reserve(nImage);
 
    for (int i = 0; i != nImage; ++i) {
        lambdaAndIndex.push_back(std::make_pair(lambda(i), i));
    }
    std::sort(lambdaAndIndex.begin(), lambdaAndIndex.end(), SortEvalueDecreasing<double>());
    //
    // Save the (sorted) eigen values
    //
    _eigenValues.clear();
    _eigenValues.reserve(nImage);
    for (int i = 0; i != nImage; ++i) {
        _eigenValues.push_back(lambdaAndIndex[i].first);
    }
    //
    // Contruct the first ncomp eigenimages in basis
    //
    int ncomp = 100;  // number of components to keep
 
    _eigenImages.clear();
    _eigenImages.reserve(ncomp < nImage ? ncomp : nImage);
 
    for (int i = 0; i < ncomp; ++i) {
        if (i >= nImage) {
            continue;
        }
 
        int const ii = lambdaAndIndex[i].second;  // the index after sorting (backwards) by eigenvalue
 
        std::shared_ptr<ImageT> eImage(new ImageT(_dimensions));
        *eImage = static_cast<typename ImageT::Pixel>(0);
 
        for (int j = 0; j != nImage; ++j) {
            int const jj = lambdaAndIndex[j].second;  // the index after sorting (backwards) by eigenvalue
            double const weight = Q(jj, ii) * (_constantWeight ? flux_bar / getFlux(jj) : 1);
            eImage->scaledPlus(weight, *_imageList[jj]);
        }
        _eigenImages.push_back(eImage);
    }
}

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 76 of file ImagePca.h.

{ return _dimensions; }

getEigenImages()

template<typename ImageT >

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

inline

Return Eigen images.

Definition at line 101 of file ImagePca.h.

{ return _eigenImages; }

getEigenValues()

template<typename ImageT >

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

inline

Return Eigen values.

Definition at line 99 of file ImagePca.h.

{ return _eigenValues; }

getImageList()

template<typename ImageT >

ImagePca< ImageT >::ImageList lsst::afw::image::ImagePca< ImageT >::getImageList

Return the list of images being analyzed.

Definition at line 86 of file ImagePca.cc.

                                                                        {
    return _imageList;
}

getMean()

template<typename ImageT >

std::shared_ptr< ImageT > lsst::afw::image::ImagePca< ImageT >::getMean

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

Definition at line 91 of file ImagePca.cc.

                                                      {
    if (_imageList.empty()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "You haven't provided any images");
    }
 
    std::shared_ptr<ImageT> mean(new ImageT(getDimensions()));
    *mean = static_cast<typename ImageT::Pixel>(0);
 
    for (typename ImageList::const_iterator ptr = _imageList.begin(), end = _imageList.end(); ptr != end;
         ++ptr) {
        *mean += **ptr;
    }
    *mean /= _imageList.size();
 
    return mean;
}

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

mask	Mask defining bad pixels
ncomp	Number 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 387 of file ImagePca.cc.

                                                                            {
    return do_updateBadPixels<ImageT>(typename ImageT::image_category(), _imageList, _fluxList, _eigenImages,
                                      mask, ncomp);
}

---

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

• /j/snowflake/release/lsstsw/stack/0.4.3/Linux64/afw/21.0.0-25-g1d57be3cd+e73869a214/include/lsst/afw/image/ImagePca.h
• /j/snowflake/release/lsstsw/stack/0.4.3/Linux64/afw/21.0.0-25-g1d57be3cd+e73869a214/src/image/ImagePca.cc