BuildSingleKernelVisitor.cc

Go to the documentation of this file.

// -*- lsst-c++ -*-
#include "boost/shared_ptr.hpp" 
#include "Eigen/Core"

#include "lsst/afw/math.h"
#include "lsst/afw/image.h"

#include "lsst/pex/exceptions/Runtime.h"
#include "lsst/pex/policy/Policy.h"
#include "lsst/pex/logging/Trace.h"

#include "lsst/ip/diffim/ImageSubtract.h"
#include "lsst/ip/diffim/KernelCandidate.h"
#include "lsst/ip/diffim/BuildSingleKernelVisitor.h"

#define DEBUG_MATRIX 0

namespace afwMath        = lsst::afw::math;
namespace afwImage       = lsst::afw::image;
namespace pexLogging     = lsst::pex::logging; 
namespace pexPolicy      = lsst::pex::policy; 
namespace pexExcept      = lsst::pex::exceptions; 
namespace ipDiffim       = lsst::ip::diffim;

namespace lsst { 
namespace ip { 
namespace diffim {
namespace detail {

    template<typename PixelT>
    BuildSingleKernelVisitor<PixelT>::BuildSingleKernelVisitor(
        lsst::afw::math::KernelList const& basisList,   
        lsst::pex::policy::Policy const& policy  
        ) :
        afwMath::CandidateVisitor(),
        _basisList(basisList),
        _policy(policy),
        _hMat(),
        _imstats(ImageStatistics<PixelT>(_policy)),
        _skipBuilt(true),
        _nRejected(0),
        _nProcessed(0),
        _useRegularization(false),
        _useCoreStats(_policy.getBool("useCoreStats")),
        _coreRadius(_policy.getInt("candidateCoreRadius"))
    {};

    template<typename PixelT>
    BuildSingleKernelVisitor<PixelT>::BuildSingleKernelVisitor(
        lsst::afw::math::KernelList const& basisList,   
        lsst::pex::policy::Policy const& policy,  
        boost::shared_ptr<Eigen::MatrixXd> hMat   
        ) :
        afwMath::CandidateVisitor(),
        _basisList(basisList),
        _policy(policy),
        _hMat(hMat),
        _imstats(ImageStatistics<PixelT>(_policy)),
        _skipBuilt(true),
        _nRejected(0),
        _nProcessed(0),
        _useRegularization(true),
        _useCoreStats(_policy.getBool("useCoreStats")),
        _coreRadius(_policy.getInt("candidateCoreRadius"))
    {};

    
    template<typename PixelT>
    void BuildSingleKernelVisitor<PixelT>::processCandidate(
        lsst::afw::math::SpatialCellCandidate *candidate
        ) {
        
        ipDiffim::KernelCandidate<PixelT> *kCandidate = 
            dynamic_cast<ipDiffim::KernelCandidate<PixelT> *>(candidate);
        if (kCandidate == NULL) {
            pexLogging::TTrace<3>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate", 
                                  "Failed to cast SpatialCellCandidate to KernelCandidate %d", 
                                  kCandidate->getId());
            throw LSST_EXCEPT(pexExcept::LogicError,
                              "Failed to cast SpatialCellCandidate to KernelCandidate");
        }
        
        if (_skipBuilt and kCandidate->isInitialized()) {
            return;
        }
        
        pexLogging::TTrace<2>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate", 
                              "Processing candidate %d", kCandidate->getId());
        pexLogging::TTrace<5>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
                              "X = %.2f Y = %.2f",
                              kCandidate->getXCenter(), 
                              kCandidate->getYCenter());
                              
        /* Build its kernel here */
        try {
            if (_useRegularization)
                kCandidate->build(_basisList, _hMat);
            else
                kCandidate->build(_basisList);

        } catch (pexExcept::Exception &e) {
            kCandidate->setStatus(afwMath::SpatialCellCandidate::BAD);
            pexLogging::TTrace<4>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate", 
                                  "Unable to process candidate %d; exception caught (%s)", 
                                  kCandidate->getId(),
                                  e.what());
            _nRejected += 1;
            return;
        } 

        if (kCandidate->getStatus() == afwMath::SpatialCellCandidate::BAD) {
            pexLogging::TTrace<4>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate", 
                                  "Candidate %d Returned BAD upon build, exiting", 
                                  kCandidate->getId());
            _nRejected += 1;
            return;
        }
            
         
         /* 
         * Make diffim and set chi2 from result.  Note that you need to use the
         * most recent kernel
         */
        MaskedImageT diffim = kCandidate->getDifferenceImage(ipDiffim::KernelCandidate<PixelT>::RECENT);
        try {
            if (_useCoreStats) 
                _imstats.apply(diffim, _coreRadius);
            else
                _imstats.apply(diffim);
        } catch (pexExcept::Exception& e) {
            pexLogging::TTrace<3>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate", 
                                  "Unable to calculate imstats for Candidate %d", kCandidate->getId()); 
            kCandidate->setStatus(afwMath::SpatialCellCandidate::BAD);
            return;
        }
        _nProcessed += 1;

        kCandidate->setChi2(_imstats.getVariance());
        
        /* When using a Pca basis, we don't reset the kernel or background,
           so we need to evaluate these locally for the Trace */
        double kSum = kCandidate->getKsum(ipDiffim::KernelCandidate<PixelT>::RECENT);
        double background = kCandidate->getBackground(ipDiffim::KernelCandidate<PixelT>::RECENT);
        
        pexLogging::TTrace<5>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate", 
                              "Chi2 = %.3f", kCandidate->getChi2());
        pexLogging::TTrace<5>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
                              "Kernel Sum = %.3f", kSum);
        pexLogging::TTrace<5>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
                              "Background = %.3f", background);
        pexLogging::TTrace<3>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
                              "Candidate %d resids = %.3f +/- %.3f sigma (%d pix)",
                              kCandidate->getId(),
                              _imstats.getMean(),
                              _imstats.getRms(),
                              _imstats.getNpix());
        
        bool meanIsNan = std::isnan(_imstats.getMean());
        bool rmsIsNan  = std::isnan(_imstats.getRms());
        if (meanIsNan || rmsIsNan) {
            kCandidate->setStatus(afwMath::SpatialCellCandidate::BAD);
            pexLogging::TTrace<4>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate", 
                                  "Rejecting candidate %d, encountered NaN",
                                  kCandidate->getId());
            _nRejected += 1;
            return;
        }
        
        if (_policy.getBool("singleKernelClipping")) {
            if (fabs(_imstats.getMean()) > _policy.getDouble("candidateResidualMeanMax")) {
                kCandidate->setStatus(afwMath::SpatialCellCandidate::BAD);
                pexLogging::TTrace<4>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate", 
                                      "Rejecting candidate %d; bad mean residual : |%.3f| > %.3f",
                                      kCandidate->getId(),
                                      _imstats.getMean(),
                                      _policy.getDouble("candidateResidualMeanMax"));
                _nRejected += 1;
            }
            else if (_imstats.getRms() > _policy.getDouble("candidateResidualStdMax")) {
                kCandidate->setStatus(afwMath::SpatialCellCandidate::BAD);
                pexLogging::TTrace<4>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate", 
                                      "Rejecting candidate %d; bad residual rms : %.3f > %.3f",
                                      kCandidate->getId(),
                                      _imstats.getRms(),
                                      _policy.getDouble("candidateResidualStdMax"));
                _nRejected += 1;
            }
            else {
                kCandidate->setStatus(afwMath::SpatialCellCandidate::GOOD);
                pexLogging::TTrace<4>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate", 
                                      "Source kernel OK");
            }
        }
        else {
            kCandidate->setStatus(afwMath::SpatialCellCandidate::GOOD);
            pexLogging::TTrace<6>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate", 
                                  "Sigma clipping not enabled");
        }
        
        /* Core resids for debugging */
        if (!(_useCoreStats)) {
            try {
                _imstats.apply(diffim, _coreRadius);
            } catch (pexExcept::Exception& e) {
                pexLogging::TTrace<3>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate", 
                                      "Unable to calculate core imstats for Candidate %d", 
                                      kCandidate->getId()); 
                kCandidate->setStatus(afwMath::SpatialCellCandidate::BAD);
                return;
            }
            pexLogging::TTrace<4>("lsst.ip.diffim.BuildSingleKernelVisitor.processCandidate",
                                  "Candidate %d core resids = %.3f +/- %.3f sigma (%d pix)",
                                  kCandidate->getId(),
                                  _imstats.getMean(),
                                  _imstats.getRms(),
                                  _imstats.getNpix());
        }
        
    }

    typedef float PixelT;

    template class BuildSingleKernelVisitor<PixelT>;

    template boost::shared_ptr<BuildSingleKernelVisitor<PixelT> >
    makeBuildSingleKernelVisitor<PixelT>(lsst::afw::math::KernelList const&,
                                         lsst::pex::policy::Policy const&);

    template boost::shared_ptr<BuildSingleKernelVisitor<PixelT> >
    makeBuildSingleKernelVisitor<PixelT>(lsst::afw::math::KernelList const&,
                                         lsst::pex::policy::Policy const&,
                                         boost::shared_ptr<Eigen::MatrixXd>);

}}}} // end of namespace lsst::ip::diffim::detail