lsst::ip::diffim::MaskedKernelSolution< InputT > Class Template Reference

#include <KernelSolution.h>

Inheritance diagram for lsst::ip::diffim::MaskedKernelSolution< InputT >:

Public Types
typedef boost::shared_ptr < MaskedKernelSolution< InputT > >	Ptr
Public Types inherited from lsst::ip::diffim::StaticKernelSolution< InputT >
typedef boost::shared_ptr < StaticKernelSolution< InputT > >	Ptr
Public Types inherited from lsst::ip::diffim::KernelSolution
enum	KernelSolvedBy {   NONE = 0, CHOLESKY_LDLT = 1, CHOLESKY_LLT = 2, LU = 3,   EIGENVECTOR = 4 }
enum	ConditionNumberType { EIGENVALUE = 0, SVD = 1 }
typedef boost::shared_ptr < KernelSolution >	Ptr
typedef lsst::afw::math::Kernel::Pixel	PixelT
typedef lsst::afw::image::Image < lsst::afw::math::Kernel::Pixel >	ImageT

Public Member Functions
	MaskedKernelSolution (lsst::afw::math::KernelList const &basisList, bool fitForBackground)
virtual	~MaskedKernelSolution ()
virtual void	buildOrig (lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate, lsst::afw::image::Mask< lsst::afw::image::MaskPixel > pixelMask)
virtual void	buildWithMask (lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate, lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &pixelMask)
virtual void	buildSingleMaskOrig (lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate, lsst::afw::geom::Box2I maskBox)
Public Member Functions inherited from lsst::ip::diffim::StaticKernelSolution< InputT >
	StaticKernelSolution (lsst::afw::math::KernelList const &basisList, bool fitForBackground)
virtual	~StaticKernelSolution ()
void	solve ()
virtual void	build (lsst::afw::image::Image< InputT > const &templateImage, lsst::afw::image::Image< InputT > const &scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &varianceEstimate)
virtual lsst::afw::math::Kernel::Ptr	getKernel ()
virtual lsst::afw::image::Image < lsst::afw::math::Kernel::Pixel > ::Ptr	makeKernelImage ()
virtual double	getBackground ()
virtual double	getKsum ()
virtual std::pair < boost::shared_ptr < lsst::afw::math::Kernel > , double >	getSolutionPair ()
Public Member Functions inherited from lsst::ip::diffim::KernelSolution
	KernelSolution (boost::shared_ptr< Eigen::MatrixXd > mMat, boost::shared_ptr< Eigen::VectorXd > bVec, bool fitForBackground)
	KernelSolution (bool fitForBackground)
	KernelSolution ()
virtual	~KernelSolution ()
virtual void	solve (Eigen::MatrixXd mMat, Eigen::VectorXd bVec)
KernelSolvedBy	getSolvedBy ()
virtual double	getConditionNumber (ConditionNumberType conditionType)
virtual double	getConditionNumber (Eigen::MatrixXd mMat, ConditionNumberType conditionType)
boost::shared_ptr < Eigen::MatrixXd >	getM ()
boost::shared_ptr < Eigen::VectorXd >	getB ()
void	printM ()
void	printB ()
void	printA ()
int	getId () const

Additional Inherited Members
Protected Member Functions inherited from lsst::ip::diffim::StaticKernelSolution< InputT >
void	_setKernel ()
	Set kernel after solution. More...
void	_setKernelUncertainty ()
	Not implemented. More...
Protected Attributes inherited from lsst::ip::diffim::StaticKernelSolution< InputT >
boost::shared_ptr < Eigen::MatrixXd >	_cMat
	K_i x R. More...
boost::shared_ptr < Eigen::VectorXd >	_iVec
	Vectorized I. More...
boost::shared_ptr < Eigen::VectorXd >	_ivVec
	Inverse variance. More...
lsst::afw::math::Kernel::Ptr	_kernel
	Derived single-object convolution kernel. More...
double	_background
	Derived differential background estimate. More...
double	_kSum
	Derived kernel sum. More...
Protected Attributes inherited from lsst::ip::diffim::KernelSolution
int	_id
	Unique ID for object. More...
boost::shared_ptr < Eigen::MatrixXd >	_mMat
	Derived least squares M matrix. More...
boost::shared_ptr < Eigen::VectorXd >	_bVec
	Derived least squares B vector. More...
boost::shared_ptr < Eigen::VectorXd >	_aVec
	Derived least squares solution matrix. More...
KernelSolvedBy	_solvedBy
	Type of algorithm used to make solution. More...
bool	_fitForBackground
	Background terms included in fit. More...
Static Protected Attributes inherited from lsst::ip::diffim::KernelSolution
static int	_SolutionId = 0
	Unique identifier for solution. More...

Detailed Description

template<typename InputT>
class lsst::ip::diffim::MaskedKernelSolution< InputT >

Definition at line 118 of file KernelSolution.h.

Member Typedef Documentation

template<typename InputT >

typedef boost::shared_ptr<MaskedKernelSolution<InputT> > lsst::ip::diffim::MaskedKernelSolution< InputT >::Ptr

Definition at line 120 of file KernelSolution.h.

Constructor & Destructor Documentation

template<typename InputT >

lsst::ip::diffim::MaskedKernelSolution< InputT >::MaskedKernelSolution	(	lsst::afw::math::KernelList const &	basisList,
		bool	fitForBackground
	)

Definition at line 503 of file KernelSolution.cc.

          :
        StaticKernelSolution<InputT>(basisList, fitForBackground)
    {};

template<typename InputT >

virtual lsst::ip::diffim::MaskedKernelSolution< InputT >::~MaskedKernelSolution ( )

inlinevirtual

Definition at line 124 of file KernelSolution.h.

{};

Member Function Documentation

template<typename InputT >

void lsst::ip::diffim::MaskedKernelSolution< InputT >::buildOrig ( lsst::afw::image::Image< InputT > const &  templateImage, lsst::afw::image::Image< InputT > const &  scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &  varianceEstimate, lsst::afw::image::Mask< lsst::afw::image::MaskPixel >  pixelMask  )

virtual

Definition at line 692 of file KernelSolution.cc.

          {

        afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
        if (varStats.getValue(afwMath::MIN) < 0.0) {
            throw LSST_EXCEPT(pexExcept::Exception, 
                              "Error: variance less than 0.0");
        }
        if (varStats.getValue(afwMath::MIN) == 0.0) {
            throw LSST_EXCEPT(pexExcept::Exception, 
                              "Error: variance equals 0.0, cannot inverse variance weight");
        }

        lsst::afw::math::KernelList basisList = 
            boost::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(this->_kernel)->getKernelList();
        std::vector<boost::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();

        /* Only BAD pixels marked in this mask */
        lsst::afw::image::Mask<lsst::afw::image::MaskPixel> sMask(pixelMask, true);
        afwImage::MaskPixel bitMask = 
            (afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("BAD") | 
             afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("SAT") |
             afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("EDGE"));
        sMask &= bitMask;
        /* TBD: Need to figure out a way to spread this; currently done in Python */

        unsigned int const nKernelParameters     = basisList.size();
        unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
        unsigned int const nParameters           = nKernelParameters + nBackgroundParameters;

        /* NOTE - we are accessing particular elements of Eigen arrays using
           these coordinates, therefore they need to be in LOCAL coordinates.
           This was written before ndarray unification.
        */
        /* Ignore known EDGE pixels for speed */
        afwGeom::Box2I shrunkLocalBBox = (*kiter)->shrinkBBox(templateImage.getBBox(afwImage::LOCAL));
        pexLog::TTrace<5>("lsst.ip.diffim.MaskedKernelSolution.build", 
                          "Limits of good pixels after convolution: %d,%d -> %d,%d (local)", 
                          shrunkLocalBBox.getMinX(), shrunkLocalBBox.getMinY(), 
                          shrunkLocalBBox.getMaxX(), shrunkLocalBBox.getMaxY());

        /* Subimages are addressed in raw pixel coordinates */
        unsigned int startCol = shrunkLocalBBox.getMinX();
        unsigned int startRow = shrunkLocalBBox.getMinY();
        unsigned int endCol   = shrunkLocalBBox.getMaxX();
        unsigned int endRow   = shrunkLocalBBox.getMaxY();
        /* Needed for Eigen block slicing operation */
        endCol += 1;
        endRow += 1;

        boost::timer t;
        t.restart();

        /* Eigen representation of input images; only the pixels that are unconvolved in cimage below */
        Eigen::MatrixXi eMask = maskToEigenMatrix(sMask).block(startRow, 
                                                               startCol, 
                                                               endRow-startRow, 
                                                               endCol-startCol);

        Eigen::MatrixXd eigenTemplate = imageToEigenMatrix(templateImage).block(startRow, 
                                                                                startCol, 
                                                                                endRow-startRow, 
                                                                                endCol-startCol);
        Eigen::MatrixXd eigenScience = imageToEigenMatrix(scienceImage).block(startRow, 
                                                                              startCol, 
                                                                              endRow-startRow, 
                                                                              endCol-startCol);
        Eigen::MatrixXd eigeniVariance = imageToEigenMatrix(varianceEstimate).block(
            startRow, startCol, endRow-startRow, endCol-startCol
        ).array().inverse().matrix();

        /* Resize into 1-D for later usage */
        eMask.resize(eMask.rows()*eMask.cols(), 1);
        eigenTemplate.resize(eigenTemplate.rows()*eigenTemplate.cols(), 1);
        eigenScience.resize(eigenScience.rows()*eigenScience.cols(), 1);
        eigeniVariance.resize(eigeniVariance.rows()*eigeniVariance.cols(), 1);

        /* Do the masking, slowly for now... */
        Eigen::MatrixXd maskedEigenTemplate(eigenTemplate.rows(), 1);
        Eigen::MatrixXd maskedEigenScience(eigenScience.rows(), 1);
        Eigen::MatrixXd maskedEigeniVariance(eigeniVariance.rows(), 1);
        int nGood = 0;
        for (int i = 0; i < eMask.rows(); i++) {
            if (eMask(i, 0) == 0) {
                maskedEigenTemplate(nGood, 0) = eigenTemplate(i, 0);
                maskedEigenScience(nGood, 0) = eigenScience(i, 0);
                maskedEigeniVariance(nGood, 0) = eigeniVariance(i, 0);
                nGood += 1;
            }
        }
        /* Can't resize() since all values are lost; use blocks */
        eigenTemplate    = maskedEigenTemplate.block(0, 0, nGood, 1);
        eigenScience     = maskedEigenScience.block(0, 0, nGood, 1);
        eigeniVariance   = maskedEigeniVariance.block(0, 0, nGood, 1);


        /* Holds image convolved with basis function */
        afwImage::Image<InputT> cimage(templateImage.getDimensions());
        
        /* Holds eigen representation of image convolved with all basis functions */
        std::vector<boost::shared_ptr<Eigen::MatrixXd> > convolvedEigenList(nKernelParameters);
        
        /* Iterators over convolved image list and basis list */
        typename std::vector<boost::shared_ptr<Eigen::MatrixXd> >::iterator eiter = 
            convolvedEigenList.begin();
        /* Create C_i in the formalism of Alard & Lupton */
        for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
            afwMath::convolve(cimage, templateImage, **kiter, false); /* cimage stores convolved image */
            
            Eigen::MatrixXd cMat = imageToEigenMatrix(cimage).block(startRow, 
                                                                    startCol, 
                                                                    endRow-startRow, 
                                                                    endCol-startCol);
            cMat.resize(cMat.rows() * cMat.cols(), 1);

            /* Do masking */
            Eigen::MatrixXd maskedcMat(cMat.rows(), 1);
            int nGood = 0;
            for (int i = 0; i < eMask.rows(); i++) {
                if (eMask(i, 0) == 0) {
                    maskedcMat(nGood, 0) = cMat(i, 0);
                    nGood += 1;
                }
            }
            cMat = maskedcMat.block(0, 0, nGood, 1);
            boost::shared_ptr<Eigen::MatrixXd> cMatPtr (new Eigen::MatrixXd(cMat));
            *eiter = cMatPtr;
        } 

        double time = t.elapsed();
        pexLog::TTrace<5>("lsst.ip.diffim.StaticKernelSolution.build", 
                          "Total compute time to do basis convolutions : %.2f s", time);
        t.restart();
        
        /* 
           Load matrix with all values from convolvedEigenList : all images
           (eigeniVariance, convolvedEigenList) must be the same size
        */
        Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
        typename std::vector<boost::shared_ptr<Eigen::MatrixXd> >::iterator eiterj = 
            convolvedEigenList.begin();
        typename std::vector<boost::shared_ptr<Eigen::MatrixXd> >::iterator eiterE = 
            convolvedEigenList.end();
        for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
            cMat.col(kidxj) = (*eiterj)->col(0);
        }
        /* Treat the last "image" as all 1's to do the background calculation. */
        if (this->_fitForBackground)
            cMat.col(nParameters-1).fill(1.);

        this->_cMat.reset(new Eigen::MatrixXd(cMat));
        this->_ivVec.reset(new Eigen::VectorXd(eigeniVariance.col(0)));
        this->_iVec.reset(new Eigen::VectorXd(eigenScience.col(0)));

        /* Make these outside of solve() so I can check condition number */
        this->_mMat.reset(
            new Eigen::MatrixXd(this->_cMat->transpose() * this->_ivVec->asDiagonal() * *(this->_cMat))
            );
        this->_bVec.reset(
            new Eigen::VectorXd(this->_cMat->transpose() * this->_ivVec->asDiagonal() * *(this->_iVec))
            );
        
    }

template<typename InputT >

void lsst::ip::diffim::MaskedKernelSolution< InputT >::buildSingleMaskOrig ( lsst::afw::image::Image< InputT > const &  templateImage, lsst::afw::image::Image< InputT > const &  scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &  varianceEstimate, lsst::afw::geom::Box2I  maskBox  )

virtual

Definition at line 863 of file KernelSolution.cc.

          {

        afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
        if (varStats.getValue(afwMath::MIN) < 0.0) {
            throw LSST_EXCEPT(pexExcept::Exception, 
                              "Error: variance less than 0.0");
        }
        if (varStats.getValue(afwMath::MIN) == 0.0) {
            throw LSST_EXCEPT(pexExcept::Exception, 
                              "Error: variance equals 0.0, cannot inverse variance weight");
        }

        lsst::afw::math::KernelList basisList = 
            boost::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(this->_kernel)->getKernelList();

        unsigned int const nKernelParameters     = basisList.size();
        unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
        unsigned int const nParameters           = nKernelParameters + nBackgroundParameters;

        std::vector<boost::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();

        /* 
           NOTE : If we are using these views in Afw's Image space, we need to
           make sure and compensate for the XY0 of the image:

           afwGeom::Box2I fullBBox = templateImage.getBBox();
           int maskStartCol = maskBox.getMinX();
           int maskStartRow = maskBox.getMinY();
           int maskEndCol   = maskBox.getMaxX();
           int maskEndRow   = maskBox.getMaxY();

        
          If we are going to be doing the slicing in Eigen matrices derived from
          the images, ignore the XY0.

           afwGeom::Box2I fullBBox = templateImage.getBBox(afwImage::LOCAL);

           int maskStartCol = maskBox.getMinX() - templateImage.getX0();
           int maskStartRow = maskBox.getMinY() - templateImage.getY0();
           int maskEndCol   = maskBox.getMaxX() - templateImage.getX0();
           int maskEndRow   = maskBox.getMaxY() - templateImage.getY0();

        */

                                                                         
        afwGeom::Box2I shrunkBBox = (*kiter)->shrinkBBox(templateImage.getBBox());

        pexLog::TTrace<5>("lsst.ip.diffim.MaskedKernelSolution.build", 
                          "Limits of good pixels after convolution: %d,%d -> %d,%d", 
                          shrunkBBox.getMinX(), shrunkBBox.getMinY(), 
                          shrunkBBox.getMaxX(), shrunkBBox.getMaxY());

        unsigned int const startCol = shrunkBBox.getMinX();
        unsigned int const startRow = shrunkBBox.getMinY();
        unsigned int const endCol   = shrunkBBox.getMaxX();
        unsigned int const endRow   = shrunkBBox.getMaxY();

        /* NOTE: no endCol/endRow += 1 for block slicing, since we are doing the
           slicing using Afw, not Eigen

           Eigen arrays have index 0,0 in the upper right, while LSST images
           have 0,0 in the lower left.  The y-axis is flipped.  When translating
           Images to Eigen matrices in ipDiffim::imageToEigenMatrix this is
           accounted for.  However, we still need to be aware of this fact if
           addressing subregions of an Eigen matrix.  This is why the slicing is
           done in Afw, its cleaner.
           
           Please see examples/maskedKernel.cc for elaboration on some of the
           tests done to make sure this indexing gets done right.

        */


        /* Inner limits; of mask box */
        int maskStartCol = maskBox.getMinX();
        int maskStartRow = maskBox.getMinY();
        int maskEndCol   = maskBox.getMaxX();
        int maskEndRow   = maskBox.getMaxY();

        /* 

        |---------------------------|
        |      Kernel Boundary      |
        |  |---------------------|  |
        |  |         Top         |  |
        |  |......_________......|  |
        |  |      |       |      |  |
        |  |  L   |  Box  |  R   |  |
        |  |      |       |      |  |
        |  |......---------......|  |
        |  |        Bottom       |  |
        |  |---------------------|  |
        |                           |
        |---------------------------|
       
        4 regions we want to extract from the pixels: top bottom left right

        */
        afwGeom::Box2I tBox = afwGeom::Box2I(afwGeom::Point2I(startCol, maskEndRow + 1),
                                             afwGeom::Point2I(endCol, endRow));
        
        afwGeom::Box2I bBox = afwGeom::Box2I(afwGeom::Point2I(startCol, startRow),
                                             afwGeom::Point2I(endCol, maskStartRow - 1));
        
        afwGeom::Box2I lBox = afwGeom::Box2I(afwGeom::Point2I(startCol, maskStartRow),
                                             afwGeom::Point2I(maskStartCol - 1, maskEndRow));
        
        afwGeom::Box2I rBox = afwGeom::Box2I(afwGeom::Point2I(maskEndCol + 1, maskStartRow),
                                             afwGeom::Point2I(endCol, maskEndRow));

        pexLog::TTrace<5>("lsst.ip.diffim.MaskedKernelSolution.build", 
                          "Upper good pixel region: %d,%d -> %d,%d", 
                          tBox.getMinX(), tBox.getMinY(), tBox.getMaxX(), tBox.getMaxY());
        pexLog::TTrace<5>("lsst.ip.diffim.MaskedKernelSolution.build", 
                          "Bottom good pixel region: %d,%d -> %d,%d", 
                          bBox.getMinX(), bBox.getMinY(), bBox.getMaxX(), bBox.getMaxY());
        pexLog::TTrace<5>("lsst.ip.diffim.MaskedKernelSolution.build", 
                          "Left good pixel region: %d,%d -> %d,%d", 
                          lBox.getMinX(), lBox.getMinY(), lBox.getMaxX(), lBox.getMaxY());
        pexLog::TTrace<5>("lsst.ip.diffim.MaskedKernelSolution.build", 
                          "Right good pixel region: %d,%d -> %d,%d", 
                          rBox.getMinX(), rBox.getMinY(), rBox.getMaxX(), rBox.getMaxY());

        std::vector<afwGeom::Box2I> boxArray;
        boxArray.push_back(tBox);
        boxArray.push_back(bBox);
        boxArray.push_back(lBox);
        boxArray.push_back(rBox);        

        int totalSize = tBox.getWidth() * tBox.getHeight();
        totalSize    += bBox.getWidth() * bBox.getHeight();
        totalSize    += lBox.getWidth() * lBox.getHeight();
        totalSize    += rBox.getWidth() * rBox.getHeight();

        Eigen::MatrixXd eigenTemplate(totalSize, 1);
        Eigen::MatrixXd eigenScience(totalSize, 1);
        Eigen::MatrixXd eigeniVariance(totalSize, 1);
        eigenTemplate.setZero();
        eigenScience.setZero();
        eigeniVariance.setZero();
        
        boost::timer t;
        t.restart();
 
        int nTerms = 0;
        typename std::vector<afwGeom::Box2I>::iterator biter = boxArray.begin();
        for (; biter != boxArray.end(); ++biter) {
            int area = (*biter).getWidth() * (*biter).getHeight();

            afwImage::Image<InputT> siTemplate(templateImage, *biter);
            afwImage::Image<InputT> siScience(scienceImage, *biter);
            afwImage::Image<InputT> sVarEstimate(varianceEstimate, *biter);

            Eigen::MatrixXd eTemplate = imageToEigenMatrix(siTemplate);
            Eigen::MatrixXd eScience = imageToEigenMatrix(siScience);
            Eigen::MatrixXd eiVarEstimate = imageToEigenMatrix(sVarEstimate).array().inverse().matrix();
            
            eTemplate.resize(area, 1);
            eScience.resize(area, 1);
            eiVarEstimate.resize(area, 1);

            eigenTemplate.block(nTerms, 0, area, 1) = eTemplate.block(0, 0, area, 1);
            eigenScience.block(nTerms, 0, area, 1) = eScience.block(0, 0, area, 1);
            eigeniVariance.block(nTerms, 0, area, 1) = eiVarEstimate.block(0, 0, area, 1);

            nTerms += area;
        }

        afwImage::Image<InputT> cimage(templateImage.getDimensions());

        std::vector<boost::shared_ptr<Eigen::MatrixXd> > convolvedEigenList(nKernelParameters);
        typename std::vector<boost::shared_ptr<Eigen::MatrixXd> >::iterator eiter = 
            convolvedEigenList.begin();
        /* Create C_i in the formalism of Alard & Lupton */
        for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
            afwMath::convolve(cimage, templateImage, **kiter, false); /* cimage stores convolved image */
            Eigen::MatrixXd cMat(totalSize, 1);
            cMat.setZero();

            int nTerms = 0;
            typename std::vector<afwGeom::Box2I>::iterator biter = boxArray.begin();
            for (; biter != boxArray.end(); ++biter) {
                int area = (*biter).getWidth() * (*biter).getHeight();

                afwImage::Image<InputT> csubimage(cimage, *biter);
                Eigen::MatrixXd esubimage = imageToEigenMatrix(csubimage);
                esubimage.resize(area, 1);
                cMat.block(nTerms, 0, area, 1) = esubimage.block(0, 0, area, 1);

                nTerms += area;
            }
            
            boost::shared_ptr<Eigen::MatrixXd> cMatPtr (new Eigen::MatrixXd(cMat));
            *eiter = cMatPtr;
            
        } 
        
        double time = t.elapsed();
        pexLog::TTrace<5>("lsst.ip.diffim.MaskedKernelSolution.build", 
                          "Total compute time to do basis convolutions : %.2f s", time);
        t.restart();

        /* 
           Load matrix with all values from convolvedEigenList : all images
           (eigeniVariance, convolvedEigenList) must be the same size
        */
        Eigen::MatrixXd cMat(eigenTemplate.col(0).size(), nParameters);
        typename std::vector<boost::shared_ptr<Eigen::MatrixXd> >::iterator eiterj = 
            convolvedEigenList.begin();
        typename std::vector<boost::shared_ptr<Eigen::MatrixXd> >::iterator eiterE = 
            convolvedEigenList.end();
        for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
            cMat.col(kidxj) = (*eiterj)->col(0);
        }
        /* Treat the last "image" as all 1's to do the background calculation. */
        if (this->_fitForBackground)
            cMat.col(nParameters-1).fill(1.);

        this->_cMat.reset(new Eigen::MatrixXd(cMat));
        this->_ivVec.reset(new Eigen::VectorXd(eigeniVariance.col(0)));
        this->_iVec.reset(new Eigen::VectorXd(eigenScience.col(0)));

        /* Make these outside of solve() so I can check condition number */
        this->_mMat.reset(
            new Eigen::MatrixXd(this->_cMat->transpose() * this->_ivVec->asDiagonal() * *(this->_cMat))
            );
        this->_bVec.reset(
            new Eigen::VectorXd(this->_cMat->transpose() * this->_ivVec->asDiagonal() * *(this->_iVec))
            );
        
    }

template<typename InputT >

void lsst::ip::diffim::MaskedKernelSolution< InputT >::buildWithMask ( lsst::afw::image::Image< InputT > const &  templateImage, lsst::afw::image::Image< InputT > const &  scienceImage, lsst::afw::image::Image< lsst::afw::image::VariancePixel > const &  varianceEstimate, lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &  pixelMask  )

virtual

Definition at line 511 of file KernelSolution.cc.

          {

        afwMath::Statistics varStats = afwMath::makeStatistics(varianceEstimate, afwMath::MIN);
        if (varStats.getValue(afwMath::MIN) < 0.0) {
            throw LSST_EXCEPT(pexExcept::Exception, 
                              "Error: variance less than 0.0");
        }
        if (varStats.getValue(afwMath::MIN) == 0.0) {
            throw LSST_EXCEPT(pexExcept::Exception, 
                              "Error: variance equals 0.0, cannot inverse variance weight");
        }

        /* Full footprint of all input images */
        afwDet::Footprint::Ptr fullFp(new afwDet::Footprint(templateImage.getBBox()));

        afwMath::KernelList basisList = 
            boost::dynamic_pointer_cast<afwMath::LinearCombinationKernel>(this->_kernel)->getKernelList();
        std::vector<boost::shared_ptr<afwMath::Kernel> >::const_iterator kiter = basisList.begin();

        /* Only BAD pixels marked in this mask */
        afwImage::MaskPixel bitMask = 
            (afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("BAD") | 
             afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("SAT") |
             afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("EDGE"));

        /* Create a Footprint that contains all the masked pixels set above */
        afwDet::Threshold threshold = afwDet::Threshold(bitMask, afwDet::Threshold::BITMASK, true);
        afwDet::FootprintSet maskFpSet(pixelMask, threshold, true);

        /* And spread it by the kernel half width */
        int growPix = (*kiter)->getCtr().getX();
        afwDet::FootprintSet maskedFpSetGrown(maskFpSet, growPix, true);
        
#if 0
        for (typename afwDet::FootprintSet::FootprintList::iterator 
                 ptr = maskedFpSetGrown.getFootprints()->begin(),
                 end = maskedFpSetGrown.getFootprints()->end(); 
             ptr != end; 
             ++ptr) {
            
            afwDet::setMaskFromFootprint(finalMask, 
                                         (**ptr),
                                         afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("BAD"));
        }
#endif

        afwImage::Mask<afwImage::MaskPixel> finalMask(pixelMask.getDimensions());
        afwDet::setMaskFromFootprintList(&finalMask, 
                                         *(maskedFpSetGrown.getFootprints()),
                                         afwImage::Mask<afwImage::MaskPixel>::getPlaneBitMask("BAD"));
        pixelMask.writeFits("pixelmask.fits");
        finalMask.writeFits("finalmask.fits");


        ndarray::Array<int, 1, 1> maskArray = 
            ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
        afwDet::flattenArray(*fullFp, finalMask.getArray(), 
                             maskArray, templateImage.getXY0()); /* Need to fake the XY0 */
        ndarray::EigenView<int, 1, 1> maskEigen(maskArray);

        ndarray::Array<InputT, 1, 1> arrayTemplate = 
            ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
        afwDet::flattenArray(*fullFp, templateImage.getArray(), 
                             arrayTemplate, templateImage.getXY0());
        ndarray::EigenView<InputT, 1, 1> eigenTemplate0(arrayTemplate);

        ndarray::Array<InputT, 1, 1> arrayScience = 
            ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
        afwDet::flattenArray(*fullFp, scienceImage.getArray(), 
                             arrayScience, scienceImage.getXY0());
        ndarray::EigenView<InputT, 1, 1> eigenScience0(arrayScience);

        ndarray::Array<afwImage::VariancePixel, 1, 1> arrayVariance = 
            ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
        afwDet::flattenArray(*fullFp, varianceEstimate.getArray(), 
                             arrayVariance, varianceEstimate.getXY0());
        ndarray::EigenView<afwImage::VariancePixel, 1, 1> eigenVariance0(arrayVariance);

        int nGood = 0;
        for (int i = 0; i < maskEigen.size(); i++) {
            if (maskEigen(i) == 0.0)
                nGood += 1;
        }

        Eigen::VectorXd eigenTemplate(nGood);
        Eigen::VectorXd eigenScience(nGood);
        Eigen::VectorXd eigenVariance(nGood);
        int nUsed = 0;
        for (int i = 0; i < maskEigen.size(); i++) {
            if (maskEigen(i) == 0.0) {
                eigenTemplate(nUsed) = eigenTemplate0(i);
                eigenScience(nUsed) = eigenScience0(i);
                eigenVariance(nUsed) = eigenVariance0(i);
                nUsed += 1;
            }
        }


        boost::timer t;
        t.restart();

        unsigned int const nKernelParameters     = basisList.size();
        unsigned int const nBackgroundParameters = this->_fitForBackground ? 1 : 0;
        unsigned int const nParameters           = nKernelParameters + nBackgroundParameters;

        /* Holds image convolved with basis function */
        afwImage::Image<InputT> cimage(templateImage.getDimensions());
        
        /* Holds eigen representation of image convolved with all basis functions */
        std::vector<boost::shared_ptr<Eigen::VectorXd> >
            convolvedEigenList(nKernelParameters);
        
        /* Iterators over convolved image list and basis list */
        typename std::vector<boost::shared_ptr<Eigen::VectorXd> >::iterator eiter = 
            convolvedEigenList.begin();

        /* Create C_i in the formalism of Alard & Lupton */
        for (kiter = basisList.begin(); kiter != basisList.end(); ++kiter, ++eiter) {
            afwMath::convolve(cimage, templateImage, **kiter, false); /* cimage stores convolved image */

            ndarray::Array<InputT, 1, 1> arrayC = 
                ndarray::allocate(ndarray::makeVector(fullFp->getArea()));
            afwDet::flattenArray(*fullFp, cimage.getArray(), 
                                 arrayC, cimage.getXY0());
            ndarray::EigenView<InputT, 1, 1> eigenC0(arrayC);

            Eigen::VectorXd eigenC(nGood);
            int nUsed = 0;
            for (int i = 0; i < maskEigen.size(); i++) {
                if (maskEigen(i) == 0.0) {
                    eigenC(nUsed) = eigenC0(i);
                    nUsed += 1;
                }
            }

            boost::shared_ptr<Eigen::VectorXd> 
                eigenCptr (new Eigen::VectorXd(eigenC));
            
            *eiter = eigenCptr;
        }
        double time = t.elapsed();
        pexLog::TTrace<5>("lsst.ip.diffim.StaticKernelSolution.buildWithMask", 
                          "Total compute time to do basis convolutions : %.2f s", time);
        t.restart();
        
        /* Load matrix with all convolved images */
        Eigen::MatrixXd cMat(eigenTemplate.size(), nParameters);
        typename std::vector<boost::shared_ptr<Eigen::VectorXd> >::iterator eiterj = 
            convolvedEigenList.begin();
        typename std::vector<boost::shared_ptr<Eigen::VectorXd> >::iterator eiterE = 
            convolvedEigenList.end();
        for (unsigned int kidxj = 0; eiterj != eiterE; eiterj++, kidxj++) {
            cMat.block(0, kidxj, eigenTemplate.size(), 1) = 
                Eigen::MatrixXd(**eiterj).block(0, 0, eigenTemplate.size(), 1);
        }
        /* Treat the last "image" as all 1's to do the background calculation. */
        if (this->_fitForBackground)
            cMat.col(nParameters-1).fill(1.);
        
        this->_cMat.reset(new Eigen::MatrixXd(cMat));
        //this->_ivVec.reset(new Eigen::VectorXd((eigenVariance.template cast<double>()).cwise().inverse()));
        //this->_iVec.reset(new Eigen::VectorXd(eigenScience.template cast<double>()));
        this->_ivVec.reset(new Eigen::VectorXd(eigenVariance.array().inverse().matrix()));
        this->_iVec.reset(new Eigen::VectorXd(eigenScience));

        /* Make these outside of solve() so I can check condition number */
        this->_mMat.reset(
            new Eigen::MatrixXd(this->_cMat->transpose() * this->_ivVec->asDiagonal() * *(this->_cMat))
            );
        this->_bVec.reset(
            new Eigen::VectorXd(this->_cMat->transpose() * this->_ivVec->asDiagonal() * *(this->_iVec))
            );
    }

---

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

• /lsst/home/lsstsw/stack/Linux64/ip_diffim/master-gcec1da163f+63/include/lsst/ip/diffim/KernelSolution.h
• /lsst/home/lsstsw/stack/Linux64/ip_diffim/master-gcec1da163f+63/src/KernelSolution.cc