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 <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 image
namespace math {
template <typename>
class MaskedVector;  // forward declaration
 
using WeightPixel = lsst::afw::image::VariancePixel;  // 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,       
    NCLIPPED = 0x8000,     
    NMASKED = 0x10000      
};
Property stringToStatisticsProperty(std::string const property);
 
class StatisticsControl {
public:
    enum WeightsBoolean { WEIGHTS_FALSE = 0, WEIGHTS_TRUE = 1, WEIGHTS_NONE };  // initial state is NONE
 
    StatisticsControl(double numSigmaClip = 3.0,  
                      int numIter = 3,            
                      lsst::afw::image::MaskPixel andMask =
                              0x0,  
                      bool isNanSafe = true,  
                      WeightsBoolean useWeights =
                              WEIGHTS_NONE  
                      )
            : _numSigmaClip(numSigmaClip),
              _numIter(numIter),
              _andMask(andMask),
              _noGoodPixelsMask(0x0),
              _isNanSafe(isNanSafe),
              _useWeights(useWeights),
              _calcErrorFromInputVariance(false),
              _calcErrorMosaicMode(false),
              _maskPropagationThresholds() {
        try {
            _noGoodPixelsMask = lsst::afw::image::Mask<>::getPlaneBitMask("NO_DATA");
        } catch (lsst::pex::exceptions::InvalidParameterError const &) {
            ;  // 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 noexcept { return _numSigmaClip; }
    int getNumIter() const noexcept { return _numIter; }
    int getAndMask() const noexcept { return _andMask; }
    int getNoGoodPixelsMask() const noexcept { return _noGoodPixelsMask; }
    bool getNanSafe() const noexcept { return _isNanSafe; }
    bool getWeighted() const noexcept { return _useWeights == WEIGHTS_TRUE ? true : false; }
    bool getWeightedIsSet() const noexcept { return _useWeights != WEIGHTS_NONE ? true : false; }
    bool getCalcErrorFromInputVariance() const noexcept { return _calcErrorFromInputVariance; }
    bool getCalcErrorMosaicMode() const noexcept { return _calcErrorMosaicMode; }
 
    void setNumSigmaClip(double numSigmaClip) {
        if (!(numSigmaClip > 0)) {
            throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                              "numSigmaClip has to be positive.");
        }
        _numSigmaClip = numSigmaClip;
    }
    void setNumIter(int numIter) {
        if (!(numIter > 0)) {
            throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                              "numIter has to be positive.");
        }
        _numIter = numIter;
    }
    void setAndMask(int andMask) { _andMask = andMask; }
    void setNoGoodPixelsMask(int noGoodPixelsMask) { _noGoodPixelsMask = noGoodPixelsMask; }
    void setNanSafe(bool isNanSafe) noexcept { _isNanSafe = isNanSafe; }
    void setWeighted(bool useWeights) noexcept { _useWeights = useWeights ? WEIGHTS_TRUE : WEIGHTS_FALSE; }
    void setCalcErrorFromInputVariance(bool calcErrorFromInputVariance) noexcept {
        _calcErrorFromInputVariance = calcErrorFromInputVariance;
    }
    void setCalcErrorMosaicMode(bool calcErrorMosaicMode) noexcept {
        _calcErrorMosaicMode = calcErrorMosaicMode;
    }
 
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
    bool _calcErrorMosaicMode;         // Calculate errors by taking mean of input variances
    std::vector<double> _maskPropagationThresholds;  // Thresholds for when to propagate mask bits,
                                                     // treated like a dict (unset bits are set to 1.0)
};
 
class Statistics final {
public:
    using Value = std::pair<double, double>;
 
    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());
 
    Statistics(Statistics const &) = default;
    Statistics(Statistics &&) = default;
    Statistics &operator=(Statistics const &) = default;
    Statistics &operator=(Statistics &&) = default;
    ~Statistics() noexcept = default;
 
    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 noexcept { 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
    int _nClipped;                                // number of pixels clipped
    int _nMasked;                                 // number of pixels masked
    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() noexcept { ; }  // never actually advance the iterator
};
template <typename ValueT>
class MaskImposter {
public:
    using x_iterator = infinite_iterator<ValueT>;
    explicit MaskImposter(ValueT val = 0) noexcept { _val[0] = val; }
    x_iterator row_begin(int) const noexcept { 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() || sctrl.getCalcErrorMosaicMode()) {
        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.getCalcErrorMosaicMode() ||
        (!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 final {
public:
    // types we'll use in Statistics
    using x_iterator = typename std::vector<ValueT>::const_iterator;
    using fast_iterator = typename std::vector<ValueT>::const_iterator;
    using Pixel = ValueT;
 
    // 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 noexcept { return _v.begin(); }
    x_iterator row_end(int) const noexcept { return _v.end(); }
    int getWidth() const noexcept { return _v.size(); }
    int getHeight() const noexcept { return 1; }
    lsst::geom::Extent2I getDimensions() const noexcept {
        return lsst::geom::Extent2I(getWidth(), getHeight());
    }
 
    bool empty() const noexcept { 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() || sctrl.getCalcErrorMosaicMode()) {
        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() || sctrl.getCalcErrorMosaicMode()) {
        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);
    }
}
}  // namespace math
}  // namespace afw
}  // namespace lsst
 
#endif