Statistics.h

Go to the documentation of this file.

// -*- LSST-C++ -*-

/*
 * LSST Data Management System
 * Copyright 2008, 2009, 2010 LSST Corporation.
 *
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the LSST License Statement and
 * the GNU General Public License along with this program.  If not,
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */

#if !defined(LSST_AFW_MATH_STATISTICS_H)
#define LSST_AFW_MATH_STATISTICS_H

#include <algorithm>
#include <cassert>
#include <limits>
#include "boost/iterator/iterator_adaptor.hpp"
#include "boost/tuple/tuple.hpp"
#include <memory>
#include "lsst/afw/image/MaskedImage.h"
#include "lsst/afw/math/MaskedVector.h"

namespace lsst {
namespace afw {
    namespace image {
        template<typename> class Image;
        template<typename, typename, typename> class MaskedImage;
    }
namespace math {
    template<typename> class MaskedVector; // forward declaration

    typedef lsst::afw::image::VariancePixel WeightPixel; // Type used for weights

enum Property {
    NOTHING = 0x0,         
    ERRORS = 0x1,          
    NPOINT = 0x2,          
    MEAN = 0x4,            
    STDEV = 0x8,           
    VARIANCE = 0x10,       
    MEDIAN = 0x20,         
    IQRANGE = 0x40,        
    MEANCLIP = 0x80,       
    STDEVCLIP = 0x100,     
    VARIANCECLIP = 0x200,  
    MIN = 0x400,           
    MAX = 0x800,           
    SUM = 0x1000,          
    MEANSQUARE = 0x2000,   
    ORMASK = 0x4000        
};
Property stringToStatisticsProperty(std::string const property);


class StatisticsControl {
    typedef enum { WEIGHTS_FALSE=0, WEIGHTS_TRUE=1, WEIGHTS_NONE } WeightsBoolean; // initial state is NONE
public:

    typedef std::shared_ptr<StatisticsControl> Ptr;
    typedef std::shared_ptr<StatisticsControl> const ConstPtr;

    StatisticsControl(
        double numSigmaClip = 3.0, 
        int numIter = 3,           
        lsst::afw::image::MaskPixel andMask = 0x0, 
        bool isNanSafe = true,     
        int useWeights = WEIGHTS_NONE      
                     ) :
        _numSigmaClip(numSigmaClip),
        _numIter(numIter),
        _andMask(andMask),
        _noGoodPixelsMask(0x0),
        _isNanSafe(isNanSafe),
        _useWeights(useWeights == 0 ? WEIGHTS_FALSE : (useWeights == 1) ? WEIGHTS_TRUE : WEIGHTS_NONE),
        _calcErrorFromInputVariance(false),
        _maskPropagationThresholds()
    {
        try {
            _noGoodPixelsMask = lsst::afw::image::Mask<>::getPlaneBitMask("NO_DATA");
        } catch(lsst::pex::exceptions::InvalidParameterError) {
            ;                           // Mask has no NO_DATA plane defined
        }

        assert(_numSigmaClip > 0);
        assert(_numIter > 0);
    }


    double getMaskPropagationThreshold(int bit) const;
    void setMaskPropagationThreshold(int bit, double threshold);

    double getNumSigmaClip() const { return _numSigmaClip; }
    int getNumIter() const { return _numIter; }
    int getAndMask() const { return _andMask; }
    int getNoGoodPixelsMask() const { return _noGoodPixelsMask; }
    bool getNanSafe() const { return _isNanSafe; }
    bool getWeighted() const { return _useWeights == WEIGHTS_TRUE ? true : false; }
    bool getWeightedIsSet() const { return _useWeights != WEIGHTS_NONE ? true : false; }
    bool getCalcErrorFromInputVariance() const { return _calcErrorFromInputVariance; }

    void setNumSigmaClip(double numSigmaClip) { assert(numSigmaClip > 0); _numSigmaClip = numSigmaClip; }
    void setNumIter(int numIter) { assert(numIter > 0); _numIter = numIter; }
    void setAndMask(int andMask) { _andMask = andMask; }
    void setNoGoodPixelsMask(int noGoodPixelsMask) { _noGoodPixelsMask = noGoodPixelsMask; }
    void setNanSafe(bool isNanSafe) { _isNanSafe = isNanSafe; }
    void setWeighted(bool useWeights) { _useWeights = useWeights ? WEIGHTS_TRUE : WEIGHTS_FALSE; }
    void setCalcErrorFromInputVariance(bool calcErrorFromInputVariance) {
        _calcErrorFromInputVariance = calcErrorFromInputVariance;
    }

private:

    friend class Statistics;

    double _numSigmaClip;                 // Number of standard deviations to clip at
    int _numIter;                         // Number of iterations
    int _andMask;                         // and-Mask to specify which mask planes to ignore
    int _noGoodPixelsMask;                // mask to set if no values are acceptable
    bool _isNanSafe;                      // Check for NaNs & Infs before running (slower)
    WeightsBoolean _useWeights;           // Calculate weighted statistics (enum because of 3-valued logic)
    bool _calcErrorFromInputVariance;     // Calculate errors from the input variances, if available
    std::vector<double> _maskPropagationThresholds; // Thresholds for when to propagate mask bits,
                                                    // treated like a dict (unset bits are set to 1.0)
};

class Statistics {
public:
    typedef std::pair<double, double> Value;

    template<typename ImageT, typename MaskT, typename VarianceT>
    explicit Statistics(ImageT const &img,
                        MaskT const &msk,
                        VarianceT const &var,
                        int const flags,
                        StatisticsControl const& sctrl = StatisticsControl());

    template<typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
    explicit Statistics(ImageT const &img,
                        MaskT const &msk,
                        VarianceT const &var,
                        WeightT const &weights,
                        int const flags,
                        StatisticsControl const& sctrl = StatisticsControl());

    Value getResult(Property const prop = NOTHING) const;

    double getError(Property const prop = NOTHING) const;
    double getValue(Property const prop = NOTHING) const;
    lsst::afw::image::MaskPixel getOrMask() const {
        return _allPixelOrMask;
    }

private:
    long _flags;                        // The desired calculation

    int _n;                             // number of pixels in the image
    Value _mean;                        // the image's mean
    Value _variance;                    // the image's variance
    double _min;                        // the image's minimum
    double _max;                        // the image's maximum
    double _sum;                        // the sum of all the image's pixels
    Value _meanclip;                    // the image's N-sigma clipped mean
    Value _varianceclip;                // the image's N-sigma clipped variance
    Value _median;                      // the image's median
    double _iqrange;                    // the image's interquartile range
    lsst::afw::image::MaskPixel _allPixelOrMask;   //  the 'or' of all masked pixels

    StatisticsControl _sctrl;           // the control structure
    bool _weightsAreMultiplicative;     // Multiply by weights rather than dividing by them

    template<typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
    void doStatistics(ImageT const &img,
                      MaskT const &msk,
                      VarianceT const &var,
                      WeightT const &weights,
                      int const flags,
                      StatisticsControl const& sctrl);
};

/*************************************  The factory functions **********************************/
template <typename ValueT>
class infinite_iterator
    : public boost::iterator_adaptor<infinite_iterator<ValueT>,
                                     const ValueT*, const ValueT,
                                     boost::forward_traversal_tag> {
public:
    infinite_iterator() : infinite_iterator::iterator_adaptor_(0) {}
    explicit infinite_iterator(const ValueT* p) : infinite_iterator::iterator_adaptor_(p) {}
private:
    friend class boost::iterator_core_access;
    void increment() { ; }              // never actually advance the iterator
};
template<typename ValueT>
class MaskImposter {
public:
    typedef infinite_iterator<ValueT> x_iterator;
    explicit MaskImposter(ValueT val = 0) { _val[0] = val; }
    x_iterator row_begin(int) const { return x_iterator(_val); }
private:
    ValueT _val[1];
};


template<typename Pixel>
Statistics makeStatistics(lsst::afw::image::Image<Pixel> const &img,
                          lsst::afw::image::Mask<image::MaskPixel> const &msk,
                          int const flags,
                          StatisticsControl const& sctrl = StatisticsControl()
                         ) {
    MaskImposter<WeightPixel> var;
    return Statistics(img, msk, var, flags, sctrl);
}


template<typename ImageT, typename MaskT, typename VarianceT>
Statistics makeStatistics(ImageT const &img,
                          MaskT const &msk,
                          VarianceT const &var,
                          int const flags,
                          StatisticsControl const& sctrl = StatisticsControl()
                         ) {
    return Statistics(img, msk, var, flags, sctrl);
}

template<typename Pixel>
Statistics makeStatistics(
        lsst::afw::image::MaskedImage<Pixel> const &mimg,
        int const flags,
        StatisticsControl const& sctrl = StatisticsControl()
                         )
{
    if (sctrl.getWeighted() || sctrl.getCalcErrorFromInputVariance()) {
        return Statistics(*mimg.getImage(), *mimg.getMask(), *mimg.getVariance(), flags, sctrl);
    } else {
        MaskImposter<WeightPixel> var;
        return Statistics(*mimg.getImage(), *mimg.getMask(), var, flags, sctrl);
    }
}

template<typename Pixel>
Statistics makeStatistics(
        lsst::afw::image::MaskedImage<Pixel> const &mimg,
        lsst::afw::image::Image<WeightPixel> const &weights,
        int const flags,
        StatisticsControl const& sctrl = StatisticsControl()
                         )
{
    if (sctrl.getWeighted() || sctrl.getCalcErrorFromInputVariance() ||
        (!sctrl.getWeightedIsSet() && (weights.getWidth() != 0 && weights.getHeight() != 0))) {
        return Statistics(*mimg.getImage(), *mimg.getMask(), *mimg.getVariance(), weights, flags, sctrl);
    } else {
        MaskImposter<WeightPixel> var;
        return Statistics(*mimg.getImage(), *mimg.getMask(), var, weights, flags, sctrl);
    }
}

Statistics makeStatistics(lsst::afw::image::Mask<lsst::afw::image::MaskPixel> const &msk,
                          int const flags,
                          StatisticsControl const& sctrl = StatisticsControl());



template<typename Pixel>
Statistics makeStatistics(
        lsst::afw::image::Image<Pixel> const &img, 
        int const flags,   
        StatisticsControl const& sctrl = StatisticsControl() 
) {
    // make a phony mask that will be compiled out
    MaskImposter<lsst::afw::image::MaskPixel> const msk;
    MaskImposter<WeightPixel> const var;
    return Statistics(img, msk, var, flags, sctrl);
}


template<typename ValueT>
class ImageImposter {
public:

    // types we'll use in Statistics
    typedef typename std::vector<ValueT>::const_iterator x_iterator;
    typedef typename std::vector<ValueT>::const_iterator fast_iterator;
    typedef ValueT Pixel;

    // constructors for std::vector<>, and copy constructor
    // These are both shallow! ... no actual copying of values
    explicit ImageImposter(std::vector<ValueT> const &v) : _v(v) { }
    explicit ImageImposter(ImageImposter<ValueT> const &img) : _v(img._getVector()) {}

    // The methods we'll use in Statistics
    x_iterator row_begin(int) const { return _v.begin(); }
    x_iterator row_end(int) const { return _v.end(); }
    int getWidth() const { return _v.size(); }
    int getHeight() const { return 1; }
    afw::geom::Extent2I getDimensions() const { return afw::geom::Extent2I(getWidth(), getHeight()); }

    bool empty() const { return _v.empty(); }
private:
    std::vector<ValueT> const &_v;                  // a private reference to the data
    std::vector<ValueT> const &_getVector() const { return _v; } // get the ref for the copyCon
};

template<typename EntryT>
Statistics makeStatistics(std::vector<EntryT> const &v, 
                          int const flags,   
                          StatisticsControl const& sctrl = StatisticsControl() 
                         ) {
    ImageImposter<EntryT> img(v);           // wrap the vector in a fake image
    MaskImposter<lsst::afw::image::MaskPixel> msk;     // instantiate a fake mask that will be compiled out.
    MaskImposter<WeightPixel> var;
    return Statistics(img, msk, var, flags, sctrl);
}

template<typename EntryT>
Statistics makeStatistics(std::vector<EntryT> const &v, 
                          std::vector<WeightPixel> const &vweights, 
                          int const flags,   
                          StatisticsControl const& sctrl = StatisticsControl() 
                         ) {
    ImageImposter<EntryT> img(v);           // wrap the vector in a fake image
    MaskImposter<lsst::afw::image::MaskPixel> msk;     // instantiate a fake mask that will be compiled out.
    MaskImposter<WeightPixel> var;

    ImageImposter<WeightPixel> weights(vweights);

    return Statistics(img, msk, var, weights, flags, sctrl);
}

template<typename EntryT>
Statistics makeStatistics(lsst::afw::math::MaskedVector<EntryT> const &mv, 
                          int const flags,   
                          StatisticsControl const& sctrl = StatisticsControl() 
                         ) {
    if (sctrl.getWeighted() || sctrl.getCalcErrorFromInputVariance()) {
        return Statistics(*mv.getImage(), *mv.getMask(), *mv.getVariance(), flags, sctrl);
    } else {
        MaskImposter<WeightPixel> var;
        return Statistics(*mv.getImage(), *mv.getMask(), var, flags, sctrl);
    }
}

template<typename EntryT>
Statistics makeStatistics(lsst::afw::math::MaskedVector<EntryT> const &mv, 
                          std::vector<WeightPixel> const &vweights, 
                          int const flags,   
                          StatisticsControl const& sctrl = StatisticsControl() 
                         ) {
    ImageImposter<WeightPixel> weights(vweights);

    if (sctrl.getWeighted() || sctrl.getCalcErrorFromInputVariance()) {
        return Statistics(*mv.getImage(), *mv.getMask(), *mv.getVariance(), weights, flags, sctrl);
    } else {
        MaskImposter<WeightPixel> var;
        return Statistics(*mv.getImage(), *mv.getMask(), var, weights, flags, sctrl);
    }
}

}}}

#endif