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

#include <ImagePca.h>

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

Public Types
typedef boost::shared_ptr< ImageT >	Ptr
typedef boost::shared_ptr < const ImageT >	ConstPtr
typedef std::vector< typename ImageT::Ptr >	ImageList

Public Member Functions
	ImagePca (bool constantWeight=true)
	ctor More...
virtual	~ImagePca ()
void	addImage (typename ImageT::Ptr img, double flux=0.0)
ImageList	getImageList () const
	Return the list of images being analyzed. More...
geom::Extent2I const	getDimensions () const
	Return the dimension of the images being analyzed. More...
ImageT::Ptr	getMean () const
virtual void	analyze ()
virtual double	updateBadPixels (unsigned long mask, int const ncomp)
std::vector< double > const &	getEigenValues () const
	Return Eigen values. More...
ImageList const &	getEigenImages () const
	Return Eigen images. More...

Private Member Functions
double	getFlux (int i) const

Private Attributes
ImageList	_imageList
std::vector< double >	_fluxList
geom::Extent2I	_dimensions
bool	_constantWeight
std::vector< double >	_eigenValues
ImageList	_eigenImages

Detailed Description

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

Definition at line 47 of file ImagePca.h.

Member Typedef Documentation

template<typename ImageT>

typedef boost::shared_ptr<const ImageT> lsst::afw::image::ImagePca< ImageT >::ConstPtr

Definition at line 50 of file ImagePca.h.

template<typename ImageT>

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

Definition at line 52 of file ImagePca.h.

template<typename ImageT>

typedef boost::shared_ptr<ImageT> lsst::afw::image::ImagePca< ImageT >::Ptr

Definition at line 49 of file ImagePca.h.

Constructor & Destructor Documentation

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 49 of file ImagePca.cc.

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

template<typename ImageT>

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

inlinevirtual

Definition at line 55 of file ImagePca.h.

{}

Member Function Documentation

template<typename ImageT>

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

Add an image to the set to be analyzed

Exceptions

lsst::pex::exceptions::LengthError

if all the images aren't the same size

Parameters

img	Image to add to set
flux	Image's flux

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);
}

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 135 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(typename ImageT::Ptr(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

        typename ImageT::Ptr 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);
    }
}

template<typename ImageT>

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

inline

Return the dimension of the images being analyzed.

Definition at line 61 of file ImagePca.h.

{ return _dimensions; }

template<typename ImageT>

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

inline

Return Eigen images.

Definition at line 70 of file ImagePca.h.

{ return _eigenImages; }

template<typename ImageT>

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

inline

Return Eigen values.

Definition at line 68 of file ImagePca.h.

{ return _eigenValues; }

template<typename ImageT>

double lsst::afw::image::ImagePca< ImageT >::getFlux ( int  i ) const

inlineprivate

Definition at line 73 of file ImagePca.h.

{ return _fluxList[i]; }

template<typename ImageT >

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

(

)

const

Return the list of images being analyzed.

Definition at line 92 of file ImagePca.cc.

                                                                        {
    return _imageList;
}

template<typename ImageT >

ImageT::Ptr lsst::afw::image::ImagePca< ImageT >::getMean

(

)

const

Return the mean of the images in ImagePca's list

Definition at line 101 of file ImagePca.cc.

                                                   {
    if (_imageList.empty()) {
        throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "You haven't provided any images");
    }

    typename ImageT::Ptr 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;
}

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

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

Parameters

mask	Mask defining bad pixels
ncomp	Number of components to use in estimate

Definition at line 417 of file ImagePca.cc.

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

Member Data Documentation

template<typename ImageT>

bool lsst::afw::image::ImagePca< ImageT >::_constantWeight

private

Definition at line 79 of file ImagePca.h.

template<typename ImageT>

geom::Extent2I lsst::afw::image::ImagePca< ImageT >::_dimensions

private

Definition at line 77 of file ImagePca.h.

template<typename ImageT>

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

private

Definition at line 82 of file ImagePca.h.

template<typename ImageT>

std::vector<double> lsst::afw::image::ImagePca< ImageT >::_eigenValues

private

Definition at line 81 of file ImagePca.h.

template<typename ImageT>

std::vector<double> lsst::afw::image::ImagePca< ImageT >::_fluxList

private

Definition at line 76 of file ImagePca.h.

template<typename ImageT>

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

private

Definition at line 75 of file ImagePca.h.

---

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

• /lsst/home/lsstsw/stack/Linux64/afw/master-gae20e2b580+3/include/lsst/afw/image/ImagePca.h
• /lsst/home/lsstsw/stack/Linux64/afw/master-gae20e2b580+3/src/image/ImagePca.cc