OpenCV中的超快速矩阵中位数（与Matlab一样快）

Question

我正在使用OpenCV编写代码，想要找到一个非常大的矩阵数组（单通道灰度，浮点数）的中值。

我尝试了几种方法，比如排序数组（使用std::sort）并选择中间条目，但与matlab中的median函数相比，速度非常慢。确切地说，matlab中需时0.25秒，在OpenCV中需要超过19秒。

我的输入图像最初是一个12位灰度图像，尺寸为3840x2748（约1050万像素），转换为float（CV_32FC1），其中所有值现在映射到范围[0,1]，在代码的某个点上，我通过调用以下函数请求中值：

double myMedianValue = medianMat(Input);

其中函数medianMat是：

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
}

我对函数medinaMat进行了时间测试，分别测试了整个函数以及各部分 - 如预期的那样，瓶颈在于：

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

这里有没有人有高效的解决方案？

谢谢！

编辑： 我尝试使用Adi Shavit回答中提供的std::nth_element。 现在，medianMat函数如下：

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];
}

运行时间从超过19秒降至3.5秒。但仍远远不及Matlab中使用中值函数的0.25秒...

- CV_User

---

试试这个算法：http://www.i-programmer.info/babbages-bag/505-quick-median.html。排序的最佳时间复杂度为*O(n log n)，但是该网站声称有一种O(n)*的中位数查找算法。 - Dan

我觉得使用OpenCV需要这么长时间很奇怪。你有考虑过将数据读入Mat所需的时间吗？你的输入和代码大小是多少？ - coincoin

你能发布你的代码吗？也许你以错误的顺序访问元素，这会导致缓存不良。 - Micka

请原谅我这个挑衅性的问题，但是找到如此大的浮点像素数组的精确中位数有什么意义呢？另一个问题：您是否了解预期值范围的任何知识？您的值跨越完整的范围+/- std::numeric_limits<float>::max()吗？提供一些关于输入数据的详细信息，也许可以将其映射到整数域，从而更容易/更快地解决此问题。 - Antonio

我已经添加了代码片段和有关我的输入图像的信息。感谢大家的评论和建议！ - CV_User

Answer 1

将数组排序并取中间元素并不是寻找中位数的最有效方法。它需要进行O(n log n)次操作。

使用C++，您应该使用std::nth_element() 函数并获取中间位置的迭代器。这是一个O(n)的操作:

nth_element 是一个 部分排序 算法，重新排列 [first, last) 中的元素，以使：

• 由指向nth的元素更改为任何该位置上出现的元素，如果 [first, last) 被排序 的话
• 在新的第n个元素之前的所有元素都小于或等于新的第n个元素之后的元素。

此外，您的原始数据是12位整数。您的实现做了一些事情，这使得与Matlab的比较具有问题:

你将数据转换为浮点格式（CV_32FC1或double或两者都有），这会消耗时间和资源。

代码中有一个额外的复制到vector<double>。

对于浮点数和尤其是双精度浮点数的操作比整数更耗费资源。

假设你的图像在内存中是连续的，OpenCV默认情况下是这样的，你应该使用CV_16C1，并直接在reshape()后的数据数组上进行操作。

另一个非常快的选择是简单地构建图像的直方图——这是对图像的单次遍历。然后，在直方图上操作，找到对应于每侧一半像素的bin——这最多是对bin的单次遍历。

OpenCV文档中有几个 教程 关于如何构建直方图。一旦你有了直方图，累加箱子的值，直到超过 3840x2748/2。这个箱子就是你的中位数。

- Adi Shavit

---

我使用了你的建议更改了 medianMat。由于评论框中的格式太有限，我将在我的回答中解决它。 - CV_User

感谢您的评论-与Matlab的比较并不是那么棘手，因为在Matlab中，完全相同的图像被转换为double类型。此外，我只计时了找到中位数的时间，并没有包括所有的转换。您能否请打出一个简短的代码片段，说明您的建议？谢谢！ - CV_User

你原来的问题使用了 std::sort，这就是我提出 nth_element 建议的原因。 - Adi Shavit

我解决了！我很快会将我的解决方案上传为答案。谢谢Adi！ - CV_User

nth_element 不是 O(n) 操作，它需要 O(n log n) 个交换。 - Roman Starkov

取决于您使用的版本。https://en.cppreference.com/w/cpp/algorithm/nth_element - Adi Shavit

Answer 2

好的。

在发布问题之前，我实际上已经尝试过这种方法，但由于一些愚蠢的错误，我将其排除为解决方案...无论如何，以下是方法：

我基本上使用4096个bin（2^12）创建了原始输入值的直方图，计算了累积分布函数（CDF）并对其进行正常化，以便它从0到1映射，并查找CDF中等于或大于0.5的最小索引。然后我将此索引除以12^2，从而找到所需的中位数值。现在它在运行时仅需要0.11秒（这是在没有重度优化的调试模式下），这比Matlab所需的时间少了一半以上。

这是函数（在我的情况下，nVals = 4096对应于12位值）:

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; }

- CV_User

---

好的，我建议的正是这个。 - Adi Shavit

@AdiShavit：回答得好。CV_User：你不需要计算所有箱子的CDF，只需计算到达0.5概率为止。 - saurabheights

Answer 3

从原始数据中找到它可能会更快。

由于原始数据只有12位值，因此只有4096种可能的不同值。这是一个很好且小的表！一次遍历所有数据，计算每个值出现的次数。这是O(n)操作。然后很容易找到中位数，只需从表的两端计数size/2项即可。

- sp2danny