KernelSumVisitor.cc

Go to the documentation of this file.

// -*- lsst-c++ -*-
#include <limits>

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

#include "lsst/ip/diffim/KernelCandidate.h"
#include "lsst/ip/diffim/KernelSumVisitor.h"

namespace afwMath        = lsst::afw::math;
namespace pexLogging     = lsst::pex::logging; 
namespace pexPolicy      = lsst::pex::policy; 
namespace pexExcept      = lsst::pex::exceptions; 

namespace lsst { 
namespace ip { 
namespace diffim {
namespace detail {
    
    template<typename PixelT>
    KernelSumVisitor<PixelT>::KernelSumVisitor(
        lsst::pex::policy::Policy const& policy 
        ) :
        afwMath::CandidateVisitor(),
        _mode(AGGREGATE),
        _kSums(std::vector<double>()),
        _kSumMean(0.),
        _kSumStd(0.),
        _dkSumMax(0.),
        _kSumNpts(0),
        _nRejected(0),
        _policy(policy) 
    {};
    
    template<typename PixelT>
    void KernelSumVisitor<PixelT>::resetKernelSum() {
        _kSums.clear();
        _kSumMean =  0.;
        _kSumStd  =  0.;
        _dkSumMax =  0.;
        _kSumNpts =  0;
        _nRejected = 0;
    }

    template<typename PixelT>
    void KernelSumVisitor<PixelT>::processCandidate(lsst::afw::math::SpatialCellCandidate 
                                                    *candidate) {

        KernelCandidate<PixelT> *kCandidate = dynamic_cast<KernelCandidate<PixelT> *>(candidate);
        if (kCandidate == NULL) {
            throw LSST_EXCEPT(pexExcept::LogicError,
                              "Failed to cast SpatialCellCandidate to KernelCandidate");
        }
        pexLogging::TTrace<6>("lsst.ip.diffim.KernelSumVisitor.processCandidate", 
                              "Processing candidate %d, mode %d", kCandidate->getId(), _mode);
        
        /* Grab all kernel sums and look for outliers */
        if (_mode == AGGREGATE) {
            _kSums.push_back(kCandidate->getKernelSolution(KernelCandidate<PixelT>::ORIG)->getKsum());
        }
        else if (_mode == REJECT) {
            if (_policy.getBool("kernelSumClipping")) {
                double kSum = 
                    kCandidate->getKernelSolution(KernelCandidate<PixelT>::ORIG)->getKsum();

                if (fabs(kSum - _kSumMean) > _dkSumMax) {
                    kCandidate->setStatus(afwMath::SpatialCellCandidate::BAD);
                    pexLogging::TTrace<4>("lsst.ip.diffim.KernelSumVisitor.processCandidate", 
                                          "Rejecting candidate %d; bad source kernel sum : (%.2f)",
                                          kCandidate->getId(),
                                          kSum);
                    _nRejected += 1;
                }
            }
            else {
                pexLogging::TTrace<6>("lsst.ip.diffim.KernelSumVisitor.processCandidate", 
                                      "Sigma clipping not enabled");
            }
        }
    }
    
    template<typename PixelT>
    void KernelSumVisitor<PixelT>::processKsumDistribution() {
        if (_kSums.size() == 0) {
            throw LSST_EXCEPT(pexExcept::Exception, 
                              "Unable to determine kernel sum; 0 candidates");
        }
        else if (_kSums.size() == 1) {
            pexLogging::TTrace<2>("lsst.ip.diffim.KernelSumVisitor.processKsumDistribution", 
                                  "WARNING: only 1 kernel candidate");
            
            _kSumMean = _kSums[0];
            _kSumStd  = 0.0;
            _kSumNpts = 1;
        }
        else {
            try {
                afwMath::Statistics stats = afwMath::makeStatistics(_kSums, 
                                                                    afwMath::NPOINT | 
                                                                    afwMath::MEANCLIP | 
                                                                    afwMath::STDEVCLIP); 
                _kSumMean = stats.getValue(afwMath::MEANCLIP);
                _kSumStd  = stats.getValue(afwMath::STDEVCLIP);
                _kSumNpts = static_cast<int>(stats.getValue(afwMath::NPOINT));
            } catch (pexExcept::Exception &e) {
                LSST_EXCEPT_ADD(e, "Unable to calculate kernel sum statistics");
                throw e;
            }
            if (std::isnan(_kSumMean)) {
                throw LSST_EXCEPT(pexExcept::Exception, 
                                  str(boost::format("Mean kernel sum returns NaN (%d points)") 
                                      % _kSumNpts));
            }
            if (std::isnan(_kSumStd)) {
                throw LSST_EXCEPT(pexExcept::Exception, 
                                  str(boost::format("Kernel sum stdev returns NaN (%d points)") 
                                      % _kSumNpts));
            }
        }
        _dkSumMax = _policy.getDouble("maxKsumSigma") * _kSumStd;
        pexLogging::TTrace<2>("lsst.ip.diffim.KernelSumVisitor.processCandidate", 
                              "Kernel Sum Distribution : %.3f +/- %.3f (%d points)", 
                              _kSumMean, _kSumStd, _kSumNpts);
    }
    
    typedef float PixelT;

    template class KernelSumVisitor<PixelT>;

    template boost::shared_ptr<KernelSumVisitor<PixelT> > 
    makeKernelSumVisitor<PixelT>(lsst::pex::policy::Policy const&);

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