Statistics.cc

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/>.
 */
 
/*
 * Support statistical operations on images
 */
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <memory>
#include <tuple>
#include <type_traits>
 
#include "lsst/pex/exceptions.h"
#include "lsst/afw/image/Image.h"
#include "lsst/afw/math/Statistics.h"
#include "lsst/geom/Angle.h"
 
using namespace std;
namespace pexExceptions = lsst::pex::exceptions;
 
namespace lsst {
namespace afw {
namespace math {
 
namespace {
double const NaN = std::numeric_limits<double>::quiet_NaN();
double const MAX_DOUBLE = std::numeric_limits<double>::max();
double const IQ_TO_STDEV = 0.741301109252802;  // 1 sigma in units of iqrange (assume Gaussian)
 
class AlwaysTrue {
public:
    template <typename T>
    bool operator()(T) const {
        return true;
    }
    template <typename Ta, typename Tb>
    bool operator()(Ta, Tb) const {
        return true;
    }
    template <typename Ta, typename Tb, typename Tc>
    bool operator()(Ta, Tb, Tc) const {
        return true;
    }
};
 
class AlwaysFalse {
public:
    template <typename T>
    bool operator()(T) const {
        return false;
    }
    template <typename Ta, typename Tb>
    bool operator()(Ta, Tb) const {
        return false;
    }
    template <typename Ta, typename Tb, typename Tc>
    bool operator()(Ta, Tb, Tc) const {
        return false;
    }
};
 
class CheckFinite {
public:
    template <typename T>
    bool operator()(T val) const {
        return std::isfinite(static_cast<float>(val));
    }
};
 
class CheckValueLtMin {
public:
    template <typename Tval, typename Tmin>
    bool operator()(Tval val, Tmin min) const {
        return (static_cast<Tmin>(val) < min);
    }
};
 
class CheckValueGtMax {
public:
    template <typename Tval, typename Tmax>
    bool operator()(Tval val, Tmax max) const {
        return (static_cast<Tmax>(val) > max);
    }
};
 
class CheckClipRange {
public:
    template <typename Tval, typename Tcen, typename Tmax>
    bool operator()(Tval val, Tcen center, Tmax cliplimit) const {
        Tmax tmp = fabs(val - center);
        return (tmp <= cliplimit);
    }
};
 
// define some abbreviated typenames for the test templates
typedef CheckFinite ChkFin;
typedef CheckValueLtMin ChkMin;
typedef CheckValueGtMax ChkMax;
typedef CheckClipRange ChkClip;
typedef AlwaysTrue AlwaysT;
typedef AlwaysFalse AlwaysF;
 
inline double varianceError(double const variance, int const n) {
    return 2 * (n - 1) * variance * variance / static_cast<double>(n * n);
}
 
typedef std::tuple<int,                // n
                   double,             // sum
                   Statistics::Value,  // mean
                   Statistics::Value,  // variance
                   double,             // min
                   double,             // max
                   image::MaskPixel    // allPixelOrMask
                   >
        StandardReturn;
 
/*
 * Functions which convert the booleans into calls to the proper templated types, one type per
 * recursion level
 */
template <typename IsFinite, typename HasValueLtMin, typename HasValueGtMax, typename InClipRange,
          bool useWeights, typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
StandardReturn processPixels(ImageT const &img, MaskT const &msk, VarianceT const &var,
                             WeightT const &weights, int const, int const nCrude, int const stride,
                             double const meanCrude, double const cliplimit,
                             bool const weightsAreMultiplicative, int const andMask,
                             bool const calcErrorFromInputVariance,
                             std::vector<double> const &maskPropagationThresholds) {
    int n = 0;
    double sumw = 0.0;   // sum(weight)  (N.b. weight will be 1.0 if !useWeights)
    double sumw2 = 0.0;  // sum(weight^2)
    double sumx = 0;     // sum(data*weight)
    double sumx2 = 0;    // sum(data*weight^2)
#if 1
    double sumvw2 = 0.0;  // sum(variance*weight^2)
#endif
    double min = (nCrude) ? meanCrude : MAX_DOUBLE;
    double max = (nCrude) ? meanCrude : -MAX_DOUBLE;
 
    image::MaskPixel allPixelOrMask = 0x0;
 
    std::vector<double> rejectedWeightsByBit(maskPropagationThresholds.size(), 0.0);
 
    for (int iY = 0; iY < img.getHeight(); iY += stride) {
        typename MaskT::x_iterator mptr = msk.row_begin(iY);
        typename VarianceT::x_iterator vptr = var.row_begin(iY);
        typename WeightT::x_iterator wptr = weights.row_begin(iY);
 
        for (typename ImageT::x_iterator ptr = img.row_begin(iY), end = ptr + img.getWidth(); ptr != end;
             ++ptr, ++mptr, ++vptr, ++wptr) {
            if (IsFinite()(*ptr) && !(*mptr & andMask) &&
                InClipRange()(*ptr, meanCrude, cliplimit)) {  // clip
 
                double const delta = (*ptr - meanCrude);
 
                if (useWeights) {
                    double weight = *wptr;
                    if (weightsAreMultiplicative) {
                        ;
                    } else {
                        if (*wptr <= 0) {
                            continue;
                        }
                        weight = 1 / weight;
                    }
 
                    sumw += weight;
                    sumw2 += weight * weight;
                    sumx += weight * delta;
                    sumx2 += weight * delta * delta;
 
                    if (calcErrorFromInputVariance) {
                        double const var = *vptr;
                        sumvw2 += var * weight * weight;
                    }
                } else {
                    sumx += delta;
                    sumx2 += delta * delta;
 
                    if (calcErrorFromInputVariance) {
                        double const var = *vptr;
                        sumvw2 += var;
                    }
                }
 
                allPixelOrMask |= *mptr;
 
                if (HasValueLtMin()(*ptr, min)) {
                    min = *ptr;
                }
                if (HasValueGtMax()(*ptr, max)) {
                    max = *ptr;
                }
                n++;
            } else {  // pixel has been clipped, rejected, etc.
                for (int bit = 0, nBits = maskPropagationThresholds.size(); bit < nBits; ++bit) {
                    image::MaskPixel mask = 1 << bit;
                    if (*mptr & mask) {
                        double weight = 1.0;
                        if (useWeights) {
                            weight = *wptr;
                            if (!weightsAreMultiplicative) {
                                if (*wptr <= 0) {
                                    continue;
                                }
                                weight = 1.0 / weight;
                            }
                        }
                        rejectedWeightsByBit[bit] += weight;
                    }
                }
            }
        }
    }
    if (n == 0) {
        min = NaN;
        max = NaN;
    }
 
    // estimate of population mean and variance.
    double mean, variance;
    if (!useWeights) {
        sumw = sumw2 = n;
    }
 
    for (int bit = 0, nBits = maskPropagationThresholds.size(); bit < nBits; ++bit) {
        double hypotheticalTotalWeight = sumw + rejectedWeightsByBit[bit];
        rejectedWeightsByBit[bit] /= hypotheticalTotalWeight;
        if (rejectedWeightsByBit[bit] > maskPropagationThresholds[bit]) {
            allPixelOrMask |= (1 << bit);
        }
    }
 
    // N.b. if sumw == 0 or sumw*sumw == sumw2 (e.g. n == 1) we'll get NaNs
    // N.b. the estimator of the variance assumes that the sample points all have the same variance;
    // otherwise, what is it that we're estimating?
    mean = sumx / sumw;
    variance = sumx2 / sumw - ::pow(mean, 2);         // biased estimator
    variance *= sumw * sumw / (sumw * sumw - sumw2);  // debias
 
    double meanVar;  // (standard error of mean)^2
    if (calcErrorFromInputVariance) {
        meanVar = sumvw2 / (sumw * sumw);
    } else {
        meanVar = variance * sumw2 / (sumw * sumw);
    }
 
    double varVar = varianceError(variance, n);  // error in variance; incorrect if useWeights is true
 
    sumx += sumw * meanCrude;
    mean += meanCrude;
 
    return StandardReturn(n, sumx, Statistics::Value(mean, meanVar), Statistics::Value(variance, varVar), min,
                          max, allPixelOrMask);
}
 
template <typename IsFinite, typename HasValueLtMin, typename HasValueGtMax, typename InClipRange,
          bool useWeights, typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
StandardReturn processPixels(ImageT const &img, MaskT const &msk, VarianceT const &var,
                             WeightT const &weights, int const flags, int const nCrude, int const stride,
                             double const meanCrude, double const cliplimit,
                             bool const weightsAreMultiplicative, int const andMask,
                             bool const calcErrorFromInputVariance, bool doGetWeighted,
                             std::vector<double> const &maskPropagationThresholds) {
    if (doGetWeighted) {
        return processPixels<IsFinite, HasValueLtMin, HasValueGtMax, InClipRange, true>(
                img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
                andMask, calcErrorFromInputVariance, maskPropagationThresholds);
    } else {
        return processPixels<IsFinite, HasValueLtMin, HasValueGtMax, InClipRange, false>(
                img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
                andMask, calcErrorFromInputVariance, maskPropagationThresholds);
    }
}
 
template <typename IsFinite, typename HasValueLtMin, typename HasValueGtMax, typename InClipRange,
          bool useWeights, typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
StandardReturn processPixels(ImageT const &img, MaskT const &msk, VarianceT const &var,
                             WeightT const &weights, int const flags, int const nCrude, int const stride,
                             double const meanCrude, double const cliplimit,
                             bool const weightsAreMultiplicative, int const andMask,
                             bool const calcErrorFromInputVariance, bool doCheckFinite, bool doGetWeighted,
                             std::vector<double> const &maskPropagationThresholds) {
    if (doCheckFinite) {
        return processPixels<CheckFinite, HasValueLtMin, HasValueGtMax, InClipRange, useWeights>(
                img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
                andMask, calcErrorFromInputVariance, doGetWeighted, maskPropagationThresholds);
    } else {
        return processPixels<AlwaysTrue, HasValueLtMin, HasValueGtMax, InClipRange, useWeights>(
                img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
                andMask, calcErrorFromInputVariance, doGetWeighted, maskPropagationThresholds);
    }
}
 
template <typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
StandardReturn getStandard(ImageT const &img, MaskT const &msk, VarianceT const &var, WeightT const &weights,
                           int const flags, bool const weightsAreMultiplicative, int const andMask,
                           bool const calcErrorFromInputVariance, bool doCheckFinite, bool doGetWeighted,
                           std::vector<double> const &maskPropagationThresholds) {
    // =====================================================
    // a crude estimate of the mean, used for numerical stability of variance
    int nCrude = 0;
    double meanCrude = 0.0;
 
    // for small numbers of values, use a small stride
    int const nPix = img.getWidth() * img.getHeight();
    int strideCrude;
    if (nPix < 100) {
        strideCrude = 2;
    } else {
        strideCrude = 10;
    }
 
    double cliplimit = -1;  // unused
    StandardReturn values = processPixels<ChkFin, AlwaysF, AlwaysF, AlwaysT, true>(
            img, msk, var, weights, flags, nCrude, strideCrude, meanCrude, cliplimit,
            weightsAreMultiplicative, andMask, calcErrorFromInputVariance, doCheckFinite, doGetWeighted,
            maskPropagationThresholds);
    nCrude = std::get<0>(values);
    double sumCrude = std::get<1>(values);
 
    meanCrude = 0.0;
    if (nCrude > 0) {
        meanCrude = sumCrude / nCrude;
    }
 
    // =======================================================
    // Estimate the full precision variance using that crude mean
    // - get the min and max as well
 
    if (flags & (MIN | MAX)) {
        return processPixels<ChkFin, ChkMin, ChkMax, AlwaysT, true>(
                img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
                andMask, calcErrorFromInputVariance, true, doGetWeighted, maskPropagationThresholds);
    } else {
        return processPixels<ChkFin, AlwaysF, AlwaysF, AlwaysT, true>(
                img, msk, var, weights, flags, nCrude, 1, meanCrude, cliplimit, weightsAreMultiplicative,
                andMask, calcErrorFromInputVariance, doCheckFinite, doGetWeighted, maskPropagationThresholds);
    }
}
 
template <typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
StandardReturn getStandard(ImageT const &img, MaskT const &msk, VarianceT const &var, WeightT const &weights,
                           int const flags, std::pair<double, double> const clipinfo,
 
                           bool const weightsAreMultiplicative, int const andMask,
                           bool const calcErrorFromInputVariance, bool doCheckFinite, bool doGetWeighted,
                           std::vector<double> const &maskPropagationThresholds) {
    double const center = clipinfo.first;
    double const cliplimit = clipinfo.second;
 
    if (std::isnan(center) || std::isnan(cliplimit)) {
        return StandardReturn(0, NaN, Statistics::Value(NaN, NaN), Statistics::Value(NaN, NaN), NaN, NaN,
                              ~0x0);
    }
 
    // =======================================================
    // Estimate the full precision variance using that crude mean
    int const stride = 1;
    int nCrude = 0;
 
    if (flags & (MIN | MAX)) {
        return processPixels<ChkFin, ChkMin, ChkMax, ChkClip, true>(
                img, msk, var, weights, flags, nCrude, stride, center, cliplimit, weightsAreMultiplicative,
                andMask, calcErrorFromInputVariance, true, doGetWeighted, maskPropagationThresholds);
    } else {  // fast loop ... just the mean & variance
        return processPixels<ChkFin, AlwaysF, AlwaysF, ChkClip, true>(
                img, msk, var, weights, flags, nCrude, stride, center, cliplimit, weightsAreMultiplicative,
                andMask, calcErrorFromInputVariance, doCheckFinite, doGetWeighted, maskPropagationThresholds);
    }
}
 
template <typename Pixel>
typename enable_if<!is_integral<Pixel>::value, double>::type 
percentile(std::vector<Pixel> &img, double const fraction) {
    assert(fraction >= 0.0 && fraction <= 1.0);
 
    int const n = img.size();
 
    if (n > 1) {
        double const idx = fraction * (n - 1);
 
        // interpolate linearly between the adjacent values
        // For efficiency:
        // - if we're asked for a fraction > 0.5,
        //    we'll do the second partial sort on shorter (upper) portion
        // - otherwise, the shorter portion will be the lower one, we'll partial-sort that.
 
        int const q1 = static_cast<int>(idx);
        int const q2 = q1 + 1;
 
        auto mid1 = img.begin() + q1;
        auto mid2 = img.begin() + q2;
        if (fraction > 0.5) {
            std::nth_element(img.begin(), mid1, img.end());
            std::nth_element(mid1, mid2, img.end());
        } else {
            std::nth_element(img.begin(), mid2, img.end());
            std::nth_element(img.begin(), mid1, mid2);
        }
 
        double val1 = static_cast<double>(*mid1);
        double val2 = static_cast<double>(*mid2);
        double w1 = (static_cast<double>(q2) - idx);
        double w2 = (idx - static_cast<double>(q1));
        return w1 * val1 + w2 * val2;
 
    } else if (n == 1) {
        return img[0];
    } else {
        return NaN;
    }
}
 
namespace {
    //
    // Helper function to estimate a floating-point quantile from integer data
    //
    // The data has been partially sorted, using nth_element
    //
    // begin:   iterator to a point which we know to be below our median
    // end:     iterator to a point which we know to be beyond our median
    // naive:   the integer value of the desired quantile
    // target:  the number of points that should be to the left of the quantile.
    //          N.b. if begin isn't the start of the data, this may not be the
    //          desired number of points.  Caveat Callor
    template<typename T>
    double
    computeQuantile(typename std::vector<T>::const_iterator begin,
                    typename std::vector<T>::const_iterator end,
                    T const naive,
                    double const target
                   )
    {
        // investigate the cumulative histogram near naive
        std::size_t left = 0;           // number of values less than naive
        std::size_t middle = 0;         // number of values equal to naive
 
        for (auto ptr = begin; ptr != end; ++ptr) {
            auto const val = *ptr;
            if (val < naive) {
                ++left;
            } else if (val == naive) {
                ++middle;
            }
        }
 
        return naive - 0.5 + (target - left)/middle;
    }
 
    template <typename Pixel>
    typename enable_if<is_integral<Pixel>::value, double>::type 
    percentile(std::vector<Pixel> &img, double const fraction) {
        assert(fraction >= 0.0 && fraction <= 1.0);
 
        auto const n = img.size();
 
        if (n == 0) {
            return NaN;
        } else if (n == 1) {
            return img[0];
        } else {
            // We need to handle ties.  The proper way to do this is to analyse the cumulative curve after
            // building the histograms (which is faster than a generic partitioning algorithm), but it's a
            // nuisance as we don't know the range of values
            //
            // This code looks clean enough, but actually the call to nth_element is expensive
            // and we *still* have to go through the array a second time
 
            double const idx = fraction*(n - 1);
 
            auto midP = img.begin() + static_cast<int>(idx);
            std::nth_element(img.begin(), midP, img.end());
            auto const naiveP = *midP;           // value of desired element
 
            return computeQuantile(img.begin(), img.end(), naiveP, fraction*n);
        }
    }
}
 
typedef std::tuple<double, double, double> MedianQuartileReturn;
namespace {
    template <typename Pixel>
    typename enable_if<!is_integral<Pixel>::value, MedianQuartileReturn>::type
    medianAndQuartiles(std::vector<Pixel> &img) {
        int const n = img.size();
 
        if (n > 1) {
            double const idx50 = 0.50 * (n - 1);
            double const idx25 = 0.25 * (n - 1);
            double const idx75 = 0.75 * (n - 1);
 
            // For efficiency:
            // - partition at 50th, then partition the two half further to get 25th and 75th
            // - to get the adjacent points (for interpolation), partition between 25/50, 50/75, 75/end
            //   these should be much smaller partitions
 
            int const q50a = static_cast<int>(idx50);
            int const q50b = q50a + 1;
            int const q25a = static_cast<int>(idx25);
            int const q25b = q25a + 1;
            int const q75a = static_cast<int>(idx75);
            int const q75b = q75a + 1;
 
            auto mid50a = img.begin() + q50a;
            auto mid50b = img.begin() + q50b;
            auto mid25a = img.begin() + q25a;
            auto mid25b = img.begin() + q25b;
            auto mid75a = img.begin() + q75a;
            auto mid75b = img.begin() + q75b;
 
            // get the 50th percentile, then get the 25th and 75th on the smaller partitions
            std::nth_element(img.begin(), mid50a, img.end());
            std::nth_element(mid50a, mid75a, img.end());
            std::nth_element(img.begin(), mid25a, mid50a);
 
            // and the adjacent points for each ... use the smallest segments available.
            std::nth_element(mid50a, mid50b, mid75a);
            std::nth_element(mid25a, mid25b, mid50a);
            std::nth_element(mid75a, mid75b, img.end());
 
            // interpolate linearly between the adjacent values
            double val50a = static_cast<double>(*mid50a);
            double val50b = static_cast<double>(*mid50b);
            double w50a = (static_cast<double>(q50b) - idx50);
            double w50b = (idx50 - static_cast<double>(q50a));
            double median = w50a * val50a + w50b * val50b;
 
            double val25a = static_cast<double>(*mid25a);
            double val25b = static_cast<double>(*mid25b);
            double w25a = (static_cast<double>(q25b) - idx25);
            double w25b = (idx25 - static_cast<double>(q25a));
            double q1 = w25a * val25a + w25b * val25b;
 
            double val75a = static_cast<double>(*mid75a);
            double val75b = static_cast<double>(*mid75b);
            double w75a = (static_cast<double>(q75b) - idx75);
            double w75b = (idx75 - static_cast<double>(q75a));
            double q3 = w75a * val75a + w75b * val75b;
 
            return MedianQuartileReturn(median, q1, q3);
        } else if (n == 1) {
            return MedianQuartileReturn(img[0], img[0], img[0]);
        } else {
            return MedianQuartileReturn(NaN, NaN, NaN);
        }
    }
 
    template <typename Pixel>
    typename enable_if<is_integral<Pixel>::value, MedianQuartileReturn>::type
    medianAndQuartiles(std::vector<Pixel> &img) {
        auto const n = img.size();
 
        if (n == 0) {
            return MedianQuartileReturn(NaN, NaN, NaN);
        } else if (n == 1) {
            return MedianQuartileReturn(img[0], img[0], img[0]);
        } else {
            // We need to handle ties.  The proper way to do this is to analyse the cumulative curve after
            // building the histograms (which is faster than a generic partitioning algorithm), but it's a
            // nuisance as we don't know the range of values
            //
            // This code looks clean enough, but actually the call to nth_element is expensive
            // and we *still* have to go through the array a second time.
        
            // For efficiency:
            // - partition at 50th, then partition the two halves further to get 25th and 75th
 
            auto mid25 = img.begin() + static_cast<int>(0.25*(n - 1));
            auto mid50 = img.begin() + static_cast<int>(0.50*(n - 1));
            auto mid75 = img.begin() + static_cast<int>(0.75*(n - 1));
 
            // get the 50th percentile, then get the 25th and 75th on the smaller partitions
            std::nth_element(img.begin(), mid50, img.end());
            std::nth_element(img.begin(), mid25, mid50);
            std::nth_element(mid50,       mid75, img.end());
 
            double const q1     = computeQuantile(img.begin(), mid50,     *mid25,
                                                  0.25*n);
            double const median = computeQuantile(mid25,       mid75,     *mid50,
                                                  0.50*n - (mid25 - img.begin()));
            double const q3     = computeQuantile(mid50,       img.end(), *mid75,
                                                  0.75*n - (mid50 - img.begin()));
 
            return MedianQuartileReturn(median, q1, q3);
        }
    }
}
 
template <typename IsFinite, typename ImageT, typename MaskT, typename VarianceT>
std::shared_ptr<std::vector<typename ImageT::Pixel> > makeVectorCopy(ImageT const &img, MaskT const &msk,
                                                                     VarianceT const &, int const andMask) {
    // Note need to keep track of allPixelOrMask here ... processPixels() does that
    // and it always gets called
    std::shared_ptr<std::vector<typename ImageT::Pixel> > imgcp(new std::vector<typename ImageT::Pixel>(0));
 
    for (int i_y = 0; i_y < img.getHeight(); ++i_y) {
        typename MaskT::x_iterator mptr = msk.row_begin(i_y);
        for (typename ImageT::x_iterator ptr = img.row_begin(i_y), end = img.row_end(i_y); ptr != end;
             ++ptr) {
            if (IsFinite()(*ptr) && !(*mptr & andMask)) {
                imgcp->push_back(*ptr);
            }
            ++mptr;
        }
    }
 
    return imgcp;
}
}  // namespace
 
double StatisticsControl::getMaskPropagationThreshold(int bit) const {
    int oldSize = _maskPropagationThresholds.size();
    if (oldSize < bit) {
        return 1.0;
    }
    return _maskPropagationThresholds[bit];
}
 
void StatisticsControl::setMaskPropagationThreshold(int bit, double threshold) {
    int oldSize = _maskPropagationThresholds.size();
    if (oldSize <= bit) {
        int newSize = bit + 1;
        _maskPropagationThresholds.resize(newSize);
        for (int i = oldSize; i < bit; ++i) {
            _maskPropagationThresholds[i] = 1.0;
        }
    }
    _maskPropagationThresholds[bit] = threshold;
}
 
Property stringToStatisticsProperty(std::string const property) {
    static std::map<std::string, Property> statisticsProperty;
    if (statisticsProperty.size() == 0) {
        statisticsProperty["NOTHING"] = NOTHING;
        statisticsProperty["ERRORS"] = ERRORS;
        statisticsProperty["NPOINT"] = NPOINT;
        statisticsProperty["MEAN"] = MEAN;
        statisticsProperty["STDEV"] = STDEV;
        statisticsProperty["VARIANCE"] = VARIANCE;
        statisticsProperty["MEDIAN"] = MEDIAN;
        statisticsProperty["IQRANGE"] = IQRANGE;
        statisticsProperty["MEANCLIP"] = MEANCLIP;
        statisticsProperty["STDEVCLIP"] = STDEVCLIP;
        statisticsProperty["VARIANCECLIP"] = VARIANCECLIP;
        statisticsProperty["MIN"] = MIN;
        statisticsProperty["MAX"] = MAX;
        statisticsProperty["SUM"] = SUM;
        statisticsProperty["MEANSQUARE"] = MEANSQUARE;
        statisticsProperty["ORMASK"] = ORMASK;
        statisticsProperty["NCLIPPED"] = NCLIPPED;
        statisticsProperty["NMASKED"] = NMASKED;
    }
    return statisticsProperty[property];
}
 
template <typename ImageT, typename MaskT, typename VarianceT>
Statistics::Statistics(ImageT const &img, MaskT const &msk, VarianceT const &var, int const flags,
                       StatisticsControl const &sctrl)
        : _flags(flags),
          _mean(NaN, NaN),
          _variance(NaN, NaN),
          _min(NaN),
          _max(NaN),
          _sum(NaN),
          _meanclip(NaN, NaN),
          _varianceclip(NaN, NaN),
          _median(NaN, NaN),
          _nClipped(0),
          _nMasked(0),
          _iqrange(NaN),
          _sctrl(sctrl),
          _weightsAreMultiplicative(false) {
    doStatistics(img, msk, var, var, _flags, _sctrl);
}
 
namespace {
template <typename T>
bool isEmpty(T const &t) {
    return t.empty();
}
 
template <typename T>
bool isEmpty(image::Image<T> const &im) {
    return (im.getWidth() == 0 && im.getHeight() == 0);
}
 
// Asserts that image dimensions are equal
template <typename ImageT1, typename ImageT2>
void checkDimensions(ImageT1 const &image1, ImageT2 const &image2) {
    if (image1.getDimensions() != image2.getDimensions()) {
        throw LSST_EXCEPT(pexExceptions::InvalidParameterError,
                          (boost::format("Image sizes don't match: %s vs %s") % image1.getDimensions() %
                           image2.getDimensions())
                                  .str());
    }
}
 
// Overloads for MaskImposter (which doesn't have a size)
template <typename ImageT, typename PixelT>
void checkDimensions(ImageT const &image1, MaskImposter<PixelT> const &image2) {}
template <typename ImageT, typename PixelT>
void checkDimensions(MaskImposter<PixelT> const &image1, ImageT const &image2) {}
}  // namespace
 
template <typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
Statistics::Statistics(ImageT const &img, MaskT const &msk, VarianceT const &var, WeightT const &weights,
                       int const flags, StatisticsControl const &sctrl)
        : _flags(flags),
          _mean(NaN, NaN),
          _variance(NaN, NaN),
          _min(NaN),
          _max(NaN),
          _sum(NaN),
          _meanclip(NaN, NaN),
          _varianceclip(NaN, NaN),
          _median(NaN, NaN),
          _nClipped(0),
          _nMasked(0),
          _iqrange(NaN),
          _sctrl(sctrl),
          _weightsAreMultiplicative(true) {
    if (!isEmpty(weights)) {
        if (_sctrl.getWeightedIsSet() && !_sctrl.getWeighted()) {
            throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                              "You must use the weights if you provide them");
        }
 
        _sctrl.setWeighted(true);
    }
    doStatistics(img, msk, var, weights, _flags, _sctrl);
}
 
template <typename ImageT, typename MaskT, typename VarianceT, typename WeightT>
void Statistics::doStatistics(ImageT const &img, MaskT const &msk, VarianceT const &var,
                              WeightT const &weights, int const flags, StatisticsControl const &sctrl) {
    int const num = img.getWidth() * img.getHeight();
    _n = num;
    if (_n == 0) {
        throw LSST_EXCEPT(pexExceptions::InvalidParameterError, "Image contains no pixels");
    }
    checkDimensions(img, msk);
    checkDimensions(img, var);
    if (sctrl.getWeighted()) {
        checkDimensions(img, weights);
    }
 
    // Check that an int's large enough to hold the number of pixels
    assert(img.getWidth() * static_cast<double>(img.getHeight()) < std::numeric_limits<int>::max());
 
    // get the standard statistics
    StandardReturn standard =
            getStandard(img, msk, var, weights, flags, _weightsAreMultiplicative, _sctrl.getAndMask(),
                        _sctrl.getCalcErrorFromInputVariance(), _sctrl.getNanSafe(), _sctrl.getWeighted(),
                        _sctrl._maskPropagationThresholds);
 
    _n = std::get<0>(standard);
    _sum = std::get<1>(standard);
    _mean = std::get<2>(standard);
    _variance = std::get<3>(standard);
    _min = std::get<4>(standard);
    _max = std::get<5>(standard);
    _allPixelOrMask = std::get<6>(standard);
 
    // ==========================================================
    // now only calculate it if it's specifically requested - these all cost more!
 
    if (flags & NMASKED) {
        _nMasked = num - _n;
    }
 
    // copy the image for any routines that will use median or quantiles
    if (flags & (MEDIAN | IQRANGE | MEANCLIP | STDEVCLIP | VARIANCECLIP)) {
        // make a vector copy of the image to get the median and quartiles (will move values)
        std::shared_ptr<std::vector<typename ImageT::Pixel> > imgcp;
        if (_sctrl.getNanSafe()) {
            imgcp = makeVectorCopy<ChkFin>(img, msk, var, _sctrl.getAndMask());
        } else {
            imgcp = makeVectorCopy<AlwaysT>(img, msk, var, _sctrl.getAndMask());
        }
 
        // if we *only* want the median, just use percentile(), otherwise use medianAndQuartiles()
        if ((flags & (MEDIAN)) && !(flags & (IQRANGE | MEANCLIP | STDEVCLIP | VARIANCECLIP))) {
            _median = Value(percentile(*imgcp, 0.5), NaN);
        } else {
            MedianQuartileReturn mq = medianAndQuartiles(*imgcp);
            _median = Value(std::get<0>(mq), NaN);
            _iqrange = std::get<2>(mq) - std::get<1>(mq);
        }
 
        if (flags & (MEANCLIP | STDEVCLIP | VARIANCECLIP)) {
            for (int i_i = 0; i_i < _sctrl.getNumIter(); ++i_i) {
                double const center = ((i_i > 0) ? _meanclip : _median).first;
                double const hwidth = (i_i > 0 && _n > 1)
                                              ? _sctrl.getNumSigmaClip() * std::sqrt(_varianceclip.first)
                                              : _sctrl.getNumSigmaClip() * IQ_TO_STDEV * _iqrange;
                std::pair<double, double> const clipinfo(center, hwidth);
 
                StandardReturn clipped = getStandard(
                        img, msk, var, weights, flags, clipinfo, _weightsAreMultiplicative,
                        _sctrl.getAndMask(), _sctrl.getCalcErrorFromInputVariance(), _sctrl.getNanSafe(),
                        _sctrl.getWeighted(), _sctrl._maskPropagationThresholds);
 
                int const nClip = std::get<0>(clipped);             // number after clipping
                _nClipped = _n - nClip;                             // number clipped
                _meanclip = std::get<2>(clipped);                   // clipped mean
                double const varClip = std::get<3>(clipped).first;  // clipped variance
 
                _varianceclip = Value(varClip, varianceError(varClip, nClip));
                // ... ignore other values
            }
        }
    }
}
 
std::pair<double, double> Statistics::getResult(Property const iProp) const {
    // if iProp == NOTHING try to return their heart's delight, as specified in the constructor
    Property const prop = static_cast<Property>(((iProp == NOTHING) ? _flags : iProp) & ~ERRORS);
 
    if (!(prop & _flags)) {  // we didn't calculate it
        throw LSST_EXCEPT(pexExceptions::InvalidParameterError,
                          (boost::format("You didn't ask me to calculate %d") % prop).str());
    }
 
    Value ret(NaN, NaN);
    switch (prop) {
        case NPOINT:
            ret.first = static_cast<double>(_n);
            if (_flags & ERRORS) {
                ret.second = 0;
            }
            break;
 
        case NCLIPPED:
            ret.first = static_cast<double>(_nClipped);
            if (_flags & ERRORS) {
                ret.second = 0;
            }
            break;
 
        case NMASKED:
            ret.first = static_cast<double>(_nMasked);
            if (_flags & ERRORS) {
                ret.second = 0;
            }
            break;
 
        case SUM:
            ret.first = static_cast<double>(_sum);
            if (_flags & ERRORS) {
                ret.second = 0;
            }
            break;
 
        // == means ==
        case MEAN:
            ret.first = _mean.first;
            if (_flags & ERRORS) {
                ret.second = ::sqrt(_mean.second);
            }
            break;
        case MEANCLIP:
            ret.first = _meanclip.first;
            if (_flags & ERRORS) {
                ret.second = ::sqrt(_meanclip.second);
            }
            break;
 
        // == stdevs & variances ==
        case VARIANCE:
            ret.first = _variance.first;
            if (_flags & ERRORS) {
                ret.second = ::sqrt(_variance.second);
            }
            break;
        case STDEV:
            ret.first = sqrt(_variance.first);
            if (_flags & ERRORS) {
                ret.second = 0.5 * ::sqrt(_variance.second) / ret.first;
            }
            break;
        case VARIANCECLIP:
            ret.first = _varianceclip.first;
            if (_flags & ERRORS) {
                ret.second = ret.second;
            }
            break;
        case STDEVCLIP:
            ret.first = sqrt(_varianceclip.first);
            if (_flags & ERRORS) {
                ret.second = 0.5 * ::sqrt(_varianceclip.second) / ret.first;
            }
            break;
 
        case MEANSQUARE:
            ret.first = (_n - 1) / static_cast<double>(_n) * _variance.first + ::pow(_mean.first, 2);
            if (_flags & ERRORS) {
                ret.second = ::sqrt(2 * ::pow(ret.first / _n, 2));  // assumes Gaussian
            }
            break;
 
        // == other stats ==
        case MIN:
            ret.first = _min;
            if (_flags & ERRORS) {
                ret.second = 0;
            }
            break;
        case MAX:
            ret.first = _max;
            if (_flags & ERRORS) {
                ret.second = 0;
            }
            break;
        case MEDIAN:
            ret.first = _median.first;
            if (_flags & ERRORS) {
                ret.second = sqrt(geom::HALFPI * _variance.first / _n);  // assumes Gaussian
            }
            break;
        case IQRANGE:
            ret.first = _iqrange;
            if (_flags & ERRORS) {
                ret.second = NaN;  // we're not estimating this properly
            }
            break;
 
        // no-op to satisfy the compiler
        case ERRORS:
            break;
        // default: redundant as 'ret' is initialized to NaN, NaN
        default:  // we must have set prop to _flags
            assert(iProp == 0);
            throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                              "getValue() may only be called without a parameter"
                              " if you asked for only one statistic");
    }
    return ret;
}
 
double Statistics::getValue(Property const prop) const { return getResult(prop).first; }
 
double Statistics::getError(Property const prop) const { return getResult(prop).second; }
 
template <>
Statistics::Statistics(image::Mask<image::MaskPixel> const &msk, image::Mask<image::MaskPixel> const &msk2,
                       image::Mask<image::MaskPixel> const &var, int const flags,
                       StatisticsControl const &sctrl)
        : _flags(flags),
          _mean(NaN, NaN),
          _variance(NaN, NaN),
          _min(NaN),
          _max(NaN),
          _meanclip(NaN, NaN),
          _varianceclip(NaN, NaN),
          _median(NaN, NaN),
          _nClipped(0),
          _iqrange(NaN),
          _sctrl(sctrl) {
    if ((flags & ~(NPOINT | SUM)) != 0x0) {
        throw LSST_EXCEPT(pexExceptions::InvalidParameterError,
                          "Statistics<Mask> only supports NPOINT and SUM");
    }
 
    typedef image::Mask<image::MaskPixel> Mask;
 
    _n = msk.getWidth() * msk.getHeight();
    if (_n == 0) {
        throw LSST_EXCEPT(pexExceptions::InvalidParameterError, "Image contains no pixels");
    }
 
    // Check that an int's large enough to hold the number of pixels
    assert(msk.getWidth() * static_cast<double>(msk.getHeight()) < std::numeric_limits<int>::max());
 
    image::MaskPixel sum = 0x0;
    for (int y = 0; y != msk.getHeight(); ++y) {
        for (Mask::x_iterator ptr = msk.row_begin(y), end = msk.row_end(y); ptr != end; ++ptr) {
            sum |= (*ptr)[0];
        }
    }
    _sum = sum;
}
 
/*
 * Although short, the definition can't be in the header as it must
 *       follow the specialization definition
 *       (g++ complained when this was in the header.)
 */
Statistics makeStatistics(image::Mask<image::MaskPixel> const &msk, int const flags,
                          StatisticsControl const &sctrl) {
    return Statistics(msk, msk, msk, flags, sctrl);
}
 
/*
 * Explicit instantiations
 *
 * explicit Statistics(MaskedImage const& img, int const flags,
 *                        StatisticsControl const& sctrl=StatisticsControl());
 */
//
#define STAT Statistics
 
typedef image::VariancePixel VPixel;
 
#define INSTANTIATE_MASKEDIMAGE_STATISTICS(TYPE)                                                       \
    template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
                              image::Image<VPixel> const &var, int const flags,                        \
                              StatisticsControl const &sctrl);                                         \
    template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
                              image::Image<VPixel> const &var, image::Image<VPixel> const &weights,    \
                              int const flags, StatisticsControl const &sctrl);                        \
    template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
                              image::Image<VPixel> const &var, ImageImposter<VPixel> const &weights,   \
                              int const flags, StatisticsControl const &sctrl)
 
#define INSTANTIATE_MASKEDIMAGE_STATISTICS_NO_MASK(TYPE)                                                \
    template STAT::Statistics(image::Image<TYPE> const &img, MaskImposter<image::MaskPixel> const &msk, \
                              image::Image<VPixel> const &var, int const flags,                         \
                              StatisticsControl const &sctrl);                                          \
    template STAT::Statistics(image::Image<TYPE> const &img, MaskImposter<image::MaskPixel> const &msk, \
                              image::Image<VPixel> const &var, image::Image<VPixel> const &weights,     \
                              int const flags, StatisticsControl const &sctrl)
 
#define INSTANTIATE_MASKEDIMAGE_STATISTICS_NO_VAR(TYPE)                                                \
    template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
                              MaskImposter<VPixel> const &var, int const flags,                        \
                              StatisticsControl const &sctrl);                                         \
    template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
                              MaskImposter<VPixel> const &var, image::Image<VPixel> const &weights,    \
                              int const flags, StatisticsControl const &sctrl);                        \
    template STAT::Statistics(image::Image<TYPE> const &img, image::Mask<image::MaskPixel> const &msk, \
                              MaskImposter<VPixel> const &var, ImageImposter<VPixel> const &weights,   \
                              int const flags, StatisticsControl const &sctrl)
 
#define INSTANTIATE_REGULARIMAGE_STATISTICS(TYPE)                                                       \
    template STAT::Statistics(image::Image<TYPE> const &img, MaskImposter<image::MaskPixel> const &msk, \
                              MaskImposter<VPixel> const &var, int const flags,                         \
                              StatisticsControl const &sctrl)
 
#define INSTANTIATE_VECTOR_STATISTICS(TYPE)                                                              \
    template STAT::Statistics(ImageImposter<TYPE> const &img, MaskImposter<image::MaskPixel> const &msk, \
                              MaskImposter<VPixel> const &var, int const flags,                          \
                              StatisticsControl const &sctrl);                                           \
    template STAT::Statistics(ImageImposter<TYPE> const &img, MaskImposter<image::MaskPixel> const &msk, \
                              MaskImposter<VPixel> const &var, ImageImposter<VPixel> const &weights,     \
                              int const flags, StatisticsControl const &sctrl)
 
#define INSTANTIATE_IMAGE_STATISTICS(T)            \
    INSTANTIATE_MASKEDIMAGE_STATISTICS(T);         \
    INSTANTIATE_MASKEDIMAGE_STATISTICS_NO_VAR(T);  \
    INSTANTIATE_MASKEDIMAGE_STATISTICS_NO_MASK(T); \
    INSTANTIATE_REGULARIMAGE_STATISTICS(T);        \
    INSTANTIATE_VECTOR_STATISTICS(T)
 
INSTANTIATE_IMAGE_STATISTICS(double);
INSTANTIATE_IMAGE_STATISTICS(float);
INSTANTIATE_IMAGE_STATISTICS(int);
INSTANTIATE_IMAGE_STATISTICS(std::uint16_t);
INSTANTIATE_IMAGE_STATISTICS(std::uint64_t);
 
}  // namespace math
}  // namespace afw
}  // namespace lsst