super fast median of matrix in OpenCV (as fast as MATLAB)

I'm writing some code in OpenCV and want to find the median value of a very large matrix array (single channel grayscale, float).

I tried several methods such as sorting the array (using std::sort()) and picking the middle entry but it is extremely slow when comparing with the median function in MATLAB. To be precise - what takes 0.25 seconds in MATLAB takes over 19 seconds in OpenCV.

My input image is originally a 12-bit greyscale image with the dimensions 3840x2748 (~10.5 megapixels), converted to float (CV_32FC1) where all the values are now mapped to the range [0,1] and at some point in the code I request the median value by calling:

double myMedianValue = medianMat(Input);

Where the function medianMat() is:

double medianMat(cv::Mat Input){    
    Input = Input.reshape(0,1); // spread Input Mat to single row
    std::vector<double> vecFromMat;
    Input.copyTo(vecFromMat); // Copy Input Mat to vector vecFromMat    
    std::sort( vecFromMat.begin(), vecFromMat.end() ); // sort vecFromMat
        if (vecFromMat.size()%2==0) {return (vecFromMat[vecFromMat.size()/2-1]+vecFromMat[vecFromMat.size()/2])/2;} // in case of even-numbered matrix
    return vecFromMat[(vecFromMat.size()-1)/2]; // odd-number of elements in matrix
}

I timed the function medianMat() by itself and also the various parts - as expected the bottleneck is in:

std::sort( vecFromMat.begin(), vecFromMat.end() ); // sort vecFromMat

Does anyone here have an efficient solution?

I have tried using std::nth_element() given in the answer of Adi Shavit.

The function medianMat() now reads as:

double medianMat(cv::Mat Input){    
    Input = Input.reshape(0,1); // spread Input Mat to single row
    std::vector<double> vecFromMat;
    Input.copyTo(vecFromMat); // Copy Input Mat to vector vecFromMat
    std::nth_element(vecFromMat.begin(), vecFromMat.begin() + vecFromMat.size() / 2, vecFromMat.end());
    return vecFromMat[vecFromMat.size() / 2];
}

The runtime has lowered from over 19 seconds to 3.5 seconds. This is still nowhere near the 0.25 second in MATLAB using the median function...

OK.

I actually tried this before posting the question and due to some silly mistakes I disqualified it as a solution... anyway here it is:

I basically create a histogram of values for my original input with 2^12 = 4096 bins, compute the CDF and normalize it so it is mapped from 0 to 1 and find the smallest index in the CDF that is equal or larger than 0.5. I then divide this index by 12^2 and thus find the median value requested. Now it runs in 0.11 seconds (and that's in debug mode without heavy optimizations) which is less than half the time required in Matlab.

Here's the function (nVals = 4096 in my case corresponding with 12-bits of values):

double medianMat(cv::Mat Input, int nVals){

// COMPUTE HISTOGRAM OF SINGLE CHANNEL MATRIX
float range[] = { 0, nVals };
const float* histRange = { range };
bool uniform = true; bool accumulate = false;
cv::Mat hist;
calcHist(&Input, 1, 0, cv::Mat(), hist, 1, &nVals, &histRange, uniform, accumulate);

// COMPUTE CUMULATIVE DISTRIBUTION FUNCTION (CDF)
cv::Mat cdf;
hist.copyTo(cdf);
for (int i = 1; i <= nVals-1; i++){
    cdf.at<float>(i) += cdf.at<float>(i - 1);
}
cdf /= Input.total();

// COMPUTE MEDIAN
double medianVal;
for (int i = 0; i <= nVals-1; i++){
    if (cdf.at<float>(i) >= 0.5) { medianVal = i;  break; }
}
return medianVal/nVals; }