BuildSingleKernelVisitor.cc

Go to the documentation of this file.

// -*- lsst-c++ -*-
 
#include <memory>
#include "Eigen/Core"
 
#include "lsst/afw/math.h"
#include "lsst/afw/image.h"
#include "lsst/log/Log.h"
#include "lsst/pex/exceptions/Runtime.h"
#include "lsst/daf/base/PropertySet.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 dafBase        = lsst::daf::base;
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::daf::base::PropertySet const& ps  
        ) :
        afwMath::CandidateVisitor(),
        _basisList(basisList),
        _ps(ps.deepCopy()),
        _hMat(),
        _imstats(ImageStatistics<PixelT>(ps)),
        _skipBuilt(true),
        _nRejected(0),
        _nProcessed(0),
        _useRegularization(false),
        _useCoreStats(ps.getAsBool("useCoreStats")),
        _coreRadius(ps.getAsInt("candidateCoreRadius"))
    {};
 
    template<typename PixelT>
    BuildSingleKernelVisitor<PixelT>::BuildSingleKernelVisitor(
        lsst::afw::math::KernelList const& basisList,   
        lsst::daf::base::PropertySet const& ps,  
        Eigen::MatrixXd const& hMat   
        ) :
        afwMath::CandidateVisitor(),
        _basisList(basisList),
        _ps(ps.deepCopy()),
        _hMat(hMat),
        _imstats(ImageStatistics<PixelT>(ps)),
        _skipBuilt(true),
        _nRejected(0),
        _nProcessed(0),
        _useRegularization(true),
        _useCoreStats(ps.getAsBool("useCoreStats")),
        _coreRadius(ps.getAsInt("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) {
            LOGL_DEBUG("TRACE2.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;
        }
 
        LOGL_DEBUG("TRACE1.ip.diffim.BuildSingleKernelVisitor.processCandidate",
                   "Processing candidate %d", kCandidate->getId());
        LOGL_DEBUG("TRACE4.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);
            LOGL_DEBUG("TRACE3.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) {
            LOGL_DEBUG("TRACE3.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) {
            LOGL_DEBUG("TRACE2.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);
 
        LOGL_DEBUG("TRACE4.ip.diffim.BuildSingleKernelVisitor.processCandidate",
                   "Chi2 = %.3f", kCandidate->getChi2());
        LOGL_DEBUG("TRACE4.ip.diffim.BuildSingleKernelVisitor.processCandidate",
                   "Kernel Sum = %.3f", kSum);
        LOGL_DEBUG("TRACE4.ip.diffim.BuildSingleKernelVisitor.processCandidate",
                   "Background = %.3f", background);
        LOGL_DEBUG("TRACE2.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);
            LOGL_DEBUG("TRACE3.ip.diffim.BuildSingleKernelVisitor.processCandidate",
                       "Rejecting candidate %d, encountered NaN",
                       kCandidate->getId());
            _nRejected += 1;
            return;
        }
 
        if (_ps->getAsBool("singleKernelClipping")) {
            if (fabs(_imstats.getMean()) > _ps->getAsDouble("candidateResidualMeanMax")) {
                kCandidate->setStatus(afwMath::SpatialCellCandidate::BAD);
                LOGL_DEBUG("TRACE3.ip.diffim.BuildSingleKernelVisitor.processCandidate",
                           "Rejecting candidate %d; bad mean residual : |%.3f| > %.3f",
                           kCandidate->getId(),
                           _imstats.getMean(),
                           _ps->getAsDouble("candidateResidualMeanMax"));
                _nRejected += 1;
            }
            else if (_imstats.getRms() > _ps->getAsDouble("candidateResidualStdMax")) {
                kCandidate->setStatus(afwMath::SpatialCellCandidate::BAD);
                LOGL_DEBUG("TRACE3.ip.diffim.BuildSingleKernelVisitor.processCandidate",
                           "Rejecting candidate %d; bad residual rms : %.3f > %.3f",
                           kCandidate->getId(),
                           _imstats.getRms(),
                           _ps->getAsDouble("candidateResidualStdMax"));
                _nRejected += 1;
            }
            else {
                kCandidate->setStatus(afwMath::SpatialCellCandidate::GOOD);
                LOGL_DEBUG("TRACE3.ip.diffim.BuildSingleKernelVisitor.processCandidate",
                           "Source kernel OK");
            }
        }
        else {
            kCandidate->setStatus(afwMath::SpatialCellCandidate::GOOD);
            LOGL_DEBUG("TRACE5.ip.diffim.BuildSingleKernelVisitor.processCandidate",
                       "Sigma clipping not enabled");
        }
 
        /* Core resids for debugging */
        if (!(_useCoreStats)) {
            try {
                _imstats.apply(diffim, _coreRadius);
            } catch (pexExcept::Exception& e) {
                LOGL_DEBUG("TRACE2.ip.diffim.BuildSingleKernelVisitor.processCandidate",
                           "Unable to calculate core imstats for Candidate %d",
                           kCandidate->getId());
                kCandidate->setStatus(afwMath::SpatialCellCandidate::BAD);
                return;
            }
            LOGL_DEBUG("TRACE3.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 std::shared_ptr<BuildSingleKernelVisitor<PixelT> >
    makeBuildSingleKernelVisitor<PixelT>(lsst::afw::math::KernelList const&,
                                         lsst::daf::base::PropertySet const&);
 
    template std::shared_ptr<BuildSingleKernelVisitor<PixelT> >
    makeBuildSingleKernelVisitor<PixelT>(lsst::afw::math::KernelList const&,
                                         lsst::daf::base::PropertySet const&,
                                         Eigen::MatrixXd const &);
 
}}}} // end of namespace lsst::ip::diffim::detail