C/C++ 中的快速 7x7 2D 中值滤波器

Question

我正在尝试将以下代码从 matlab 转换为 c++

function data = process(data)
    data = medfilt2(data, [7 7], 'symmetric');
    mask = fspecial('gaussian', [35 35], 12);
    data = imfilter(data, mask, 'replicate', 'same');
    maximum = max(data(:));
    data = 1 ./ ( data/maximum );
    data(data > 10) = 16;
end

我在 medfilt2 中的问题 - 这是一个 2D 中值滤波器，我需要它来支持每像素 10 位和更多图像。

1.我研究过openCV，它有一个支持16位的5x5中值滤波器，但7x7只支持字节。

http://docs.opencv.org/2.4/modules/imgproc/doc/filtering.html?highlight=medianblur#medianblur

2.我也看过英特尔 IPP，但我只能看到一维中值滤波器 https://software.intel.com/en-us/node/502283

二维过滤器是否有快速实现？
寻找类似的东西：

1. Fast Median Search: An ANSI C Implementation 使用并行编程和矢量化 (AVX/SSE) 操作...
2. 二维数字信号处理 II。变换和中值滤波器。 由 T.S.Huang 编辑。施普林格出版社。 1981 年。

FAST MEDIAN FILTERING WITH IMPLEMENTATIONS IN C/C++/C#/VB.NET/Delphi中有更多代码示例。

我还找到了Median Filtering in Constant Time。

Answer 1

由于 OpenCV 没有为大内核大小（大于 5）实现 16 位中值滤波器这一事实，我尝试了三种不同的策略。

都是基于Huang的[2]滑动窗口算法。也就是说，当窗口从左向右滑动时，通过删除和插入像素条目来更新直方图。这对于 8 位图像非常简单，并且已经在 OpenCV 中实现。但是，65536 bin 的直方图比较大，计算起来有点困难。

...算法仍然保持 O(log r)，但存储考虑使其对于 16 位图像不切实际，对于浮点图像也不可能。 [3]

我在适用的情况下使用了algorithm C++ 标准库，并没有实现 Weiss 的额外优化策略。

1) 一个简单的排序实现。我认为这是任意像素类型（尤其是浮点数）的最佳起点。

// copy pixels in the sliding window to a temporary vec and
// compute the median value (size is always odd)
memcpy( &v[0], &window[0], window.size() * sizeof(_Type) );
std::vector< _Type >::iterator it = v.begin() + v.size()/2;
std::nth_element( v.begin(), it, v.end() );
return *it;

2) 稀疏直方图。我们不想跨过 65536 个 bin 来找到每个像素的中值，那么如何存储稀疏直方图呢？同样，这适用于所有像素类型，但如果窗口中的所有像素都不同（例如浮动），则没有意义。

typedef std::map< _Type, int > Map;
//...
// inside the sliding window, update the histogram as follows
for ( /* pixels to remove */ )
{
    // _Type px
    Map::iterator it = map.find( px );
    if ( it->second > 1 )
        it->second -= 1;
    else
        map.erase( it );
}
// ...
for ( /* pixels to add */ )
{
    // _Type px
    Map::iterator lower = map.lower_bound( px );
    if ( lower != map.end() && lower->first == px )
        lower->second += 1;
    else
        map.insert( lower, std::pair<_Type,int>( px, 1 ) );
}
//... and compute the median by integrating from the one end until
// until the appropriate sum is reached ..

3) 密集的直方图。所以这是密集的直方图，但不是简单的 65536 数组，而是通过将其划分为子 bin 来简化搜索，例如：

[0...65535] <- px
[0...4095] <- px / 16
[0...255] <- px / 256
[0...15] <- px / 4096

这会使插入速度变慢（按恒定时间），但搜索速度要快得多。我发现 16 是个不错的数字。

图我测试了方法 (1) red、(2) blue 和 (3) black 相互对比以及 8bpp OpenCV（绿色）。除了 OpenCV，输入图像都是 16-bpp 灰度。虚线在动态范围 [0,255] 处被截断，平滑线在 [0, 8020] 处被截断（通过乘以 16 和平滑以增加像素值的更多方差）。

有趣的是稀疏直方图随着像素值方差的增加而发散。第 N 个元素始终是安全的赌注，OpenCV 是最快的（如果 8bpp 可以的话），而密集的直方图则落后。

我使用的是 Windows 7、8 x 3.4 GHz 和 Visual Studio v. 10。我的是运行多线程的，OpenCV 实现是单线程的。输入图片大小为 2136x3201（http://i.imgur.com/gg9Z2aB.jpg，来自 Vogue）。

[2]：Huang，T：“二维信号处理II：变换 和中值滤波器”，1981

[3]：Weiss，B：“快速中值和双边滤波”，2006 年

Answer 2

我在网上找到了这个，它和 OpenCV 的算法一样，但是扩展到 16 位并针对 SSE 进行了优化。

https://github.com/aoles/EBImage/blob/master/src/medianFilter.c