Find average w/o extreme values

Question

I know there is an option to find the average of a Mat in OpenCV:

cv::mean(mat);

I want to know if in OpenCV there is also an option to find the average without the extreme values (e.g. just the values between 10% to 90%).

Answer 1

I dont know OpenCV, but I doubt that it has a ready to use function for this. However, a naive implementation could look like this:

double m = cv::mean(mat);
Mat temp = mat;
... set all elements in temp to 0, where abs(temp[i][j] - m) > tolerance
... and count those elements in count
int N = mat.total(); // total number of elements
m = cv::sum(temp) / (N-count)

EDIT: Actually this is not exactly what the question was asking for. However, if one can assume a gaussian distribution of the values, one could estimate the value of tolerance based on the standard deviation (has to be computed) to exclude the upper/lower 10% of the data.

Answer 2

No, there is no OpenCV function to do this. You can, however, implement your own.

---

The trickiest part is to compute the values that correspond to your percentages. This can be easily achieved computing the cumulative histogram of the image.

However, to make the approach general, you can't know which values are in the matrix, so you need rely on maps.

Note that if you are working only on CV_8U images, you can optimize knowing that you'll have at most 256 different values. To implement this, you can have an hint here.

So this is a possible implementation, that works on Mat with at most 4 channels (as cv::mean), and without knowing a priori the possible number of different values. You can check commenting / decommenting parts in the example matrix intialization that it performs correctly:

#include <opencv2/opencv.hpp>
#include <vector>
#include <numeric>
#include <map>
#include <iostream>
using namespace cv;
using namespace std;

double robustMeanC1(const Mat1d& src, Vec2d bounds)
{
    // Compute histogram (with no predefined range)
    map<double, int> hist;
    for (int r = 0; r < src.rows; ++r)
    {
        for (int c = 0; c < src.cols; ++c)
        {
            double key = src(r,c);
            if (hist.count(key) == 0) {
                hist.insert(make_pair(key, 1));
            }
            else {
                hist[key]++;
            }
        }
    }

    // Get vectors from map
    vector<double> vals;
    vector<int> sums;
    vals.reserve(hist.size());
    sums.reserve(hist.size());
    for (auto kv : hist)
    {
        vals.push_back(kv.first);
        sums.push_back(kv.second);
    }

    // Compute cumulative histogram
    vector<int> cumhist(sums);
    for (int i=1; i<sums.size(); ++i)
    {
        cumhist[i] = cumhist[i - 1] + sums[i];
    }

    // Compute bounds
    int total = src.rows * src.cols;
    double low_bound = (total * bounds[0]) / 100.0;
    double upp_bound = (total * bounds[1]) / 100.0;
    int low_index = distance(cumhist.begin(), upper_bound(cumhist.begin(), cumhist.end(), low_bound));
    int upp_index = distance(cumhist.begin(), upper_bound(cumhist.begin(), cumhist.end(), upp_bound));

    if (low_index == upp_index) {upp_index++;}

    // Compute mean
    double mean = 0.0;
    int count = 0;
    for (int i = low_index; i < upp_index; ++i)
    {
        mean += vals[i] * sums[i];
        count += sums[i];
    }
    mean /= count;

    return mean;
}

Scalar robustMean(const Mat& src, Vec2d bounds) 
{
    Mat m;
    src.convertTo(m, CV_64F);

    Scalar res(0.0, 0.0, 0.0, 0.0);

    if (m.channels() == 1)
    {
        res[0] = robustMeanC1(m, bounds);
    } 
    else
    {
        vector<Mat1d> planes;
        split(m, planes);

        if (planes.size() > 4)
        {
            // Error, at most 4 channels
            return Scalar(0,0,0,0);
        }

        for (int i = 0; i < planes.size(); ++i)
        {
            res[i] = robustMeanC1(planes[i], bounds);
        }
    }
    return res;
}



int main()
{
    Mat1d m(10,10, 5.f);
    m(Range(0,1), Range::all()) = 2.0;
    //m(Range(1, 2), Range::all()) = 80.0;
    //randu(m, Scalar(0), Scalar(1));

    //Mat3b m = imread("path_to_image");

    Scalar rs = robustMean(m, Vec2d(10, 90));
    Scalar s = mean(m);

    cout << "Robust Mean: " << rs << endl;
    cout << "       Mean: " << s << endl;

    return 0;
}

Answer 3

I would simply sort the Mat elements and take the mean of the truncated vector

#include <algorithm>
#include <vector>
// c in [0,1] the portion of middvalues added to the mean
template<class _T> _T avg( std::vector< _T > & vec, double c )
{
    if ( c < 0.0 )
        c = 0.0;
    else if ( c > 1.0 )
        c = 1.0;
    const size_t len = (size_t)( c * (double)vec.size() );
    if ( len == 0 )
        return 0.0;
    std::vector< _T >::iterator beg = vec.begin();
    std::vector< _T >::iterator end = vec.end();
    if ( len < vec.size() )
    {    
        beg += ( vec.size() - len )/2;
        end = beg + len;
        std::nth_element( vec.begin(), beg, vec.end() );
        std::nth_element( beg, end, vec.end() );
    }
    double sum = 0.0, d = 0.0;
    for ( std::vector<_T>::iterator it = beg; it!=end; ++it, d+=1.0 )
        sum += *it;
    return sum/d;
}
// fill the vector and compute for each channel separately.

Here tail and head the same portion for simplicity.