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

我正在尝试将以下代码从 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 https://en.wikipedia.org/wiki/OpenCV它有一个 5x5 中值滤波器，支持 16 位，但 7x7 只支持字节。

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

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

是否有快速实现 2D 滤波器的方法？

我正在寻找类似的东西：

1. 快速中值搜索：ANSI C 实现 http://ndevilla.free.fr/median使用并行编程和矢量化（AVX https://en.wikipedia.org/wiki/Advanced_Vector_Extensions/SSE https://en.wikipedia.org/wiki/Streaming_SIMD_Extensions） 运营...
2. 二维数字信号处理 II。变换和中值滤波器。 黄 T.S. 编辑。施普林格出版社。 1981年。

还有更多代码示例在 C/C++/C#/VB.NET/Delphi 中实现的快速中值滤波 http://www.sergejusz.com/engineering_tips/median_filter.htm.

我还发现恒定时间中值滤波 http://nomis80.org/ctmf.html.

由于 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 数组，而是通过将其划分为子箱来使搜索变得更容易，例如：

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

这使得插入速度稍慢（恒定时间），但搜索速度却快得多。我发现 16 是一个很好的数字。

Figure我测试了方法 (1) 红色、(2) 蓝色和 (3) 黑色相互对比以及 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 (https://i.stack.imgur.com/CsGpg.jpg https://i.stack.imgur.com/CsGpg.jpg，来自《时尚》）。

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

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