OpenCV 中的滤波函数
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了OpenCV 中的滤波函数相关的知识，希望对你有一定的参考价值。
参考技术A也称为 box filter、均值滤波器，就是简单地将每个像素的值替换成邻域平均值。
如果用 kernel （也称为 mask) 表示，就是
如果采用积分图的方法，可以更快的计算这种 box filter 的结果。
在积分图中，只需要三次加法运算，一次乘法运算即可，即通过积分图，算出 kernel 内部区域的像素和，然后取平均。
积分图中每一点 的值是原图中对应位置左上角区域所有值的和：
积分图的计算可以很高效：
每次计算只需要新增一个像素值，其他值都是之前已经计算出来的。
积分图一旦计算出来，对任意矩形区域内 像素和 的计算都可以在常数时间（即计算时间固定，与区域的大小无关）内完成，例如：
在高斯滤波器中，当前像素值取邻域的 加权平均 ，离当前像素越近，权重越大，权重服从高斯分布。
在实际应用中，几乎总是首选高斯滤波器，很少直接用 box filter.
上述命令中，最后两个参数 kernel size 和 如果只设置一个，则另一个可以通过以下公式推出：
第二个式子很好理解，就是借助高斯函数的性质（距离均值 3 个标准差范围内的取值占总数的 99.7%），因此窗口大小就是 3 倍的 *2 （两边）然后再加上 1 （自身）。
第一个式子与第二个非常近似，但是又做了一些微调。
上述高斯滤波器内部实际上是先调用如下函数，产生服从高斯分布的一系列权重：
上述 9 个权重是经过归一化的，即和为 1，其公式为
其中 满足归一化的要求，ksize 必须是奇数。如果要生成二维的高斯矩阵权重，则是先产生一个权重列向量，然后令该列向量与自身的转置相乘，得到高斯矩阵权重，最后再统一进行归一化，保证矩阵所有元素和为 1。
另外，还可以分别产生两个不同的高斯权重向量，对应行和列方向上的高斯模糊权重，然后把它们相乘得到高斯矩阵。由于满足这种分离的性质，高斯滤波器被称为可分离的滤波器。
前边在进行滤波操作时，都只包含线性操作（算数平均、加权平均）。
另外还可以采用非线性操作，对应非线性滤波器。非线性滤波器不能表示成 kernel 矩阵卷积操作的形式。
中值滤波器是一种非线性滤波器。它把当前像素和邻域像素组成一个集合，排序之后，选择中间值（即排序中间位置的数值）替换当前像素值。
椒盐噪声 ：像素随机替换成白色或者黑点。在通讯时，如果部分像素值丢失，就会产生这种噪声。
中值滤波器可以有效的消除椒盐噪声，因为这些噪声点在排序时很难成为中间值，因此全都被剔除了。
Sobel 也是线性滤波器，只不过采用了特殊的 kernel 矩阵：
分别针对水平方向和垂直方向的操作。
用上述 kernel 进行操作，就是计算水平或者垂直方向像素值的差分，近似反映了像素值水平和垂直变化的速度，合在一起就是图像梯度的近似。
在默认情况下，差分运算的结果很可能超过 [0,255] 这个范围，而且有正有负，应该用 CV_16S 数据类型表示。经过上述缩放和偏移之后，才勉强适合用 CV_8U 表示，但还是需要饱和截断操作。
在分别得到横向、纵向变化率之后，可以整合起来计算梯度的大小
一般如果要显示最后的 sobel 边缘检测的结果，还需要把上述模值转化到 [0,255] 范围内。
sobel 实际上包含了平滑和求导两个操作，其中邻域像素累加相当于高斯平滑，距离越近的像素权重越高。
sobel 的 kernel size 可以选择 1, 3, 5 和 7，其中 1 代表 1×3 或者 3×1，此时是没有高斯平滑的。
对于大的 size，这种平滑更明显。此时，sobel 不是高通滤波器，而是带通滤波器，既消除了部分高频，又消除了部分低频。
与 Sobel 算子类似的还有其他几个计算梯度的算子，只是采用不同的 kernel.
上述所有的滤波器都是近似计算图像函数的一阶导数，像素变化大的区域计算得到的值较大，平坦的区域计算值较小。
sobel 算子通过对图片函数求导，那些数值绝对值较高的点对应了边界区域：
如果继续求二阶导，则导数较大的点对应了过零点：
因此，也可以通过搜索二阶导的过零点来检测边界点。
Laplacian 算子的定义 ：
对照 Hessian 矩阵：
Laplacian 算子实际上就是 Hessian 矩阵的 Trace。
具体到图像操作中，二阶导有如下表达式：
所以最终 Laplacian 算子表达式为：
在具体实现中，可以用以下 kernel 进行卷积操作：
Laplacian 算子具有旋转 90 度不变性，即一幅图旋转 90 度及其倍数，对应的 Laplacian 算子操作结果相同。
为了得到更好的旋转不变性，可以将 Laplacian 算子 kernel 扩展为如下形式：
这样就具有了旋转 45 度及其倍数的不变性。
Laplacian 算子对噪声比较敏感，因此一般在进行 Laplacian 之前先进行高斯平滑滤波。
两个步骤合并称为 LoG (Laplacian of Gaussian)。
在具体实现中，我们并不需要先高斯再拉普拉斯，而是两步并作一步：将拉普拉斯算子作用在高斯 kernel 上，得到新的 kernel，再与 image 做卷积：
最后作用在 位置上的卷积权重为
同样也是通过 设定滤波范围。
对高斯函数取拉普拉斯算子操作是什么样子的？
二维情况下得到的曲面很像“墨西哥草帽”。
的大小决定了检测的粗粒度：
Difference of Gaussians
为了减少 LoG 计算量，用两个不同 的高斯做差，来近似 LoG
上图中两个 的取值好像反了。。。
C++ 和 OpenCV 中的 SSE 均值滤波器
【中文标题】C++ 和 OpenCV 中的 SSE 均值滤波器【英文标题】：SSE mean filter in c++ and OpenCV 【发布时间】：2021-05-15 14:11:52 【问题描述】：我想修改 OpenCV 均值滤波器的代码以使用 Intel 内在函数。我是 SSE 新手，我真的不知道从哪里开始。我在网上查了很多资源，但都没有成功。
这是程序：
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
using namespace cv;
using namespace std;
int main()
int A[3][3] = 1, 1, 1 , 1, 1, 1 , 1, 1, 1 ;
int c = 0;
int d = 0;
Mat var1 = imread("images.jpg", 1);
Mat var2(var1.rows, var1.cols, CV_8UC3, Scalar(0, 0, 0));
for (int i = 0; i < var1.rows; i++)
var2.at<Vec3b>(i, 0) = var1.at<Vec3b>(i, 0);
var2.at<Vec3b>(i, var1.cols - 1) = var1.at<Vec3b>(i, var1.cols - 1);
for (int i = 0; i < var1.cols; i++)
var2.at<Vec3b>(0, i) = var1.at<Vec3b>(0, i);
var2.at<Vec3b>(var1.rows - 1, i) = var1.at<Vec3b>(var1.rows - 1, i);
for (int i = 0; i < var1.rows; i++)
for (int j = 0; j < var1.cols; j++)
c = 0;
for (int m = i; m < var1.rows; m++, c++)
if (c < 3)
d = 0;
for (int n = j; n < var1.cols; n++, d++)
if (d < 3)
if ((i + 1) < var1.rows && (j + 1) < var1.cols)
var2.at<Vec3b>(i + 1, j + 1)[0] += var1.at<Vec3b>(m, n)[0] * A[m - i][n - j] / 9;
var2.at<Vec3b>(i + 1, j + 1)[1] += var1.at<Vec3b>(m, n)[1] * A[m - i][n - j] / 9;
var2.at<Vec3b>(i + 1, j + 1)[2] += var1.at<Vec3b>(m, n)[2] * A[m - i][n - j] / 9;
imshow("window1", var1);
imshow("window2", var2);
waitKey(0);
return(0);
我发现困难的部分是理解如何转换最里面的 2 个循环，其中计算平均值。任何帮助将不胜感激。
【问题讨论】：
你首先需要从一个更高效的标量实现开始——上面的代码效率非常低，向量化一个低效的实现是没有意义的。先去掉所有的冗余，再看SIMD。 您能指出您认为上面的代码效率低下的地方吗？ 两个大问题是：(a) 将 9 个像素乘以 1，因此对于 RGB，每个像素有 27 次冗余乘法运算，以及 (b) 对重叠像素进行冗余求和 - 对于您可以维持的平均内核部分和（例如列），并且每个新像素仅对一列求和。您应该能够实现一个 3x3 均值滤波器，每个通道每个像素有 4 个加法和一个乘法（定点除以 9）。 你试过blur()吗？
代码不好。在您的内部循环中，您是否迭代直到 var1.cols 仅获得 3 个值。更好的方法是 const int count = std::min(3, var1.rows - i)
【参考方案1】：
只是为了好玩，我认为从 3x3 均值滤波器的简单实现开始，然后逐步优化，最终实现 SIMD (SSE)，测量每个阶段的吞吐量改进可能会很有趣。
1 - Mean_3_3_ref - 参考实现
这只是一个简单的标量实现，我们将使用它作为吞吐量和验证进一步实现的基准：
void Mean_3_3_ref(const Mat &image_in, Mat &image_out)
for (int y = 1; y < image_in.rows - 1; ++y)
for (int x = 1; x < image_in.cols - 1; ++x)
for (int c = 0; c < 3; ++c)
image_out.at<Vec3b>(y, x)[c] = (image_in.at<Vec3b>(y - 1, x - 1)[c] +
image_in.at<Vec3b>(y - 1, x )[c] +
image_in.at<Vec3b>(y - 1, x + 1)[c] +
image_in.at<Vec3b>(y , x - 1)[c] +
image_in.at<Vec3b>(y , x )[c] +
image_in.at<Vec3b>(y , x + 1)[c] +
image_in.at<Vec3b>(y + 1, x - 1)[c] +
image_in.at<Vec3b>(y + 1, x )[c] +
image_in.at<Vec3b>(y + 1, x + 1)[c] + 4) / 9;
2 - Mean_3_3_scalar - 一些优化的标量实现
利用对连续列求和的冗余 - 我们保存最后两列的总和，以便我们只需要在每次迭代中计算一个新的列总和（每个通道）：
void Mean_3_3_scalar(const Mat &image_in, Mat &image_out)
for (int y = 1; y < image_in.rows - 1; ++y)
int r_1, g_1, b_1;
int r0, g0, b0;
int r1, g1, b1;
r_1 = g_1 = b_1 = 0;
r0 = g0 = b0 = 0;
for (int yy = y - 1; yy <= y + 1; ++yy)
r_1 += image_in.at<Vec3b>(yy, 0)[0];
g_1 += image_in.at<Vec3b>(yy, 0)[1];
b_1 += image_in.at<Vec3b>(yy, 0)[2];
r0 += image_in.at<Vec3b>(yy, 1)[0];
g0 += image_in.at<Vec3b>(yy, 1)[1];
b0 += image_in.at<Vec3b>(yy, 1)[2];
for (int x = 1; x < image_in.cols - 1; ++x)
r1 = g1 = b1 = 0;
for (int yy = y - 1; yy <= y + 1; ++yy)
r1 += image_in.at<Vec3b>(yy, x + 1)[0];
g1 += image_in.at<Vec3b>(yy, x + 1)[1];
b1 += image_in.at<Vec3b>(yy, x + 1)[2];
image_out.at<Vec3b>(y, x)[0] = (r_1 + r0 + r1 + 4) / 9;
image_out.at<Vec3b>(y, x)[1] = (g_1 + g0 + g1 + 4) / 9;
image_out.at<Vec3b>(y, x)[2] = (b_1 + b0 + b1 + 4) / 9;
r_1 = r0;
g_1 = g0;
b_1 = b0;
r0 = r1;
g0 = g1;
b0 = b1;
3 - Mean_3_3_scalar_opt - 进一步优化的标量实现
根据Mean_3_3_scalar，还可以通过缓存指向我们正在处理的每一行的指针来消除 OpenCV 开销：
void Mean_3_3_scalar_opt(const Mat &image_in, Mat &image_out)
for (int y = 1; y < image_in.rows - 1; ++y)
const uint8_t * const input_1 = image_in.ptr(y - 1);
const uint8_t * const input0 = image_in.ptr(y);
const uint8_t * const input1 = image_in.ptr(y + 1);
uint8_t * const output = image_out.ptr(y);
int r_1 = input_1[0] + input0[0] + input1[0];
int g_1 = input_1[1] + input0[1] + input1[1];
int b_1 = input_1[2] + input0[2] + input1[2];
int r0 = input_1[3] + input0[3] + input1[3];
int g0 = input_1[4] + input0[4] + input1[4];
int b0 = input_1[5] + input0[5] + input1[5];
for (int x = 1; x < image_in.cols - 1; ++x)
int r1 = input_1[x * 3 + 3] + input0[x * 3 + 3] + input1[x * 3 + 3];
int g1 = input_1[x * 3 + 4] + input0[x * 3 + 4] + input1[x * 3 + 4];
int b1 = input_1[x * 3 + 5] + input0[x * 3 + 5] + input1[x * 3 + 5];
output[x * 3 ] = (r_1 + r0 + r1 + 4) / 9;
output[x * 3 + 1] = (g_1 + g0 + g1 + 4) / 9;
output[x * 3 + 2] = (b_1 + b0 + b1 + 4) / 9;
r_1 = r0;
g_1 = g0;
b_1 = b0;
r0 = r1;
g0 = g1;
b0 = b1;
4 - Mean_3_3_blur - 利用 OpenCV 的模糊功能
OpenCV 有一个名为blur 的函数，它基于函数boxFilter，这只是均值滤波器的另一个名称。由于多年来 OpenCV 代码已经进行了相当大的优化（在许多情况下使用 SIMD），让我们看看这是否对我们的标量代码有很大的改进：
void Mean_3_3_blur(const Mat &image_in, Mat &image_out)
blur(image_in, image_out, Size(3, 3));
5 - Mean_3_3_SSE - SSE 实施
这是一个相当有效的 SIMD 实现。它使用与上述标量代码相同的技术来消除处理连续像素时的冗余：
#include <tmmintrin.h> // Note: requires SSSE3 (aka MNI)
inline void Load2(const ssize_t offset, const uint8_t* const src, __m128i& vh, __m128i& vl)
const __m128i v = _mm_loadu_si128((__m128i *)(src + offset));
vh = _mm_unpacklo_epi8(v, _mm_setzero_si128());
vl = _mm_unpackhi_epi8(v, _mm_setzero_si128());
inline void Store2(const ssize_t offset, uint8_t* const dest, const __m128i vh, const __m128i vl)
__m128i v = _mm_packus_epi16(vh, vl);
_mm_storeu_si128((__m128i *)(dest + offset), v);
template <int SHIFT> __m128i ShiftL(const __m128i v0, const __m128i v1) return _mm_alignr_epi8(v1, v0, SHIFT * sizeof(short));
template <int SHIFT> __m128i ShiftR(const __m128i v0, const __m128i v1) return _mm_alignr_epi8(v1, v0, 16 - SHIFT * sizeof(short));
template <int CHANNELS> void Mean_3_3_SSE_Impl(const Mat &image_in, Mat &image_out)
const int nx = image_in.cols;
const int ny = image_in.rows;
const int kx = 3 / 2; // x, y borders
const int ky = 3 / 2;
const int kScale = 3 * 3; // scale factor = total number of pixels in sum
const __m128i vkScale = _mm_set1_epi16((32768 + kScale / 2) / kScale);
const int nx0 = ((nx + kx) * CHANNELS + 15) & ~15; // round up total width to multiple of 16
int x, y;
for (y = ky; y < ny - ky; ++y)
const uint8_t * const input_1 = image_in.ptr(y - 1);
const uint8_t * const input0 = image_in.ptr(y);
const uint8_t * const input1 = image_in.ptr(y + 1);
uint8_t * const output = image_out.ptr(y);
__m128i vsuml_1, vsumh0, vsuml0;
__m128i vh, vl;
vsuml_1 = _mm_set1_epi16(0);
Load2(0, input_1, vsumh0, vsuml0);
Load2(0, input0, vh, vl);
vsumh0 = _mm_add_epi16(vsumh0, vh);
vsuml0 = _mm_add_epi16(vsuml0, vl);
Load2(0, input1, vh, vl);
vsumh0 = _mm_add_epi16(vsumh0, vh);
vsuml0 = _mm_add_epi16(vsuml0, vl);
for (x = 0; x < nx0; x += 16)
__m128i vsumh1, vsuml1, vsumh, vsuml;
Load2((x + 16), input_1, vsumh1, vsuml1);
Load2((x + 16), input0, vh, vl);
vsumh1 = _mm_add_epi16(vsumh1, vh);
vsuml1 = _mm_add_epi16(vsuml1, vl);
Load2((x + 16), input1, vh, vl);
vsumh1 = _mm_add_epi16(vsumh1, vh);
vsuml1 = _mm_add_epi16(vsuml1, vl);
vsumh = _mm_add_epi16(vsumh0, ShiftR<CHANNELS>(vsuml_1, vsumh0));
vsuml = _mm_add_epi16(vsuml0, ShiftR<CHANNELS>(vsumh0, vsuml0));
vsumh = _mm_add_epi16(vsumh, ShiftL<CHANNELS>(vsumh0, vsuml0));
vsuml = _mm_add_epi16(vsuml, ShiftL<CHANNELS>(vsuml0, vsumh1));
// round mean
vsumh = _mm_mulhrs_epi16(vsumh, vkScale);
vsuml = _mm_mulhrs_epi16(vsuml, vkScale);
Store2(x, output, vsumh, vsuml);
vsuml_1 = vsuml0;
vsumh0 = vsumh1;
vsuml0 = vsuml1;
void Mean_3_3_SSE(const Mat &image_in, Mat &image_out)
const int channels = image_in.channels();
switch (channels)
case 1:
Mean_3_3_SSE_Impl<1>(image_in, image_out);
break;
case 3:
Mean_3_3_SSE_Impl<3>(image_in, image_out);
break;
default:
throw("Unsupported format.");
break;
结果
我在 2.4 GHz 的第 8 代 Core i9 (MacBook Pro 16,1) 上对上述所有实现进行了基准测试，图像大小为 2337 行 x 3180 列。编译器是 Apple clang 版本 12.0.5 (clang-1205.0.22.9)，唯一的优化开关是 -O3。 OpenCV 版本是 4.5.0（通过Homebrew）。 （注意：我验证了 Mean_3_3_blur 的 cv::blur 函数已分派到 AVX2 实现。）结果：
Mean_3_3_ref 62153 µs
Mean_3_3_scalar 41144 µs = 1.51062x
Mean_3_3_scalar_opt 26238 µs = 2.36882x
Mean_3_3_blur 20121 µs = 3.08896x
Mean_3_3_SSE 4838 µs = 12.84680x
注意事项
我在所有实现中都忽略了边框像素 - 如果需要，可以使用原始图像中的像素或使用其他形式的边缘像素处理来填充这些像素。
代码不是“工业实力”——它只是为了基准测试目的而编写的。
还有一些可能的优化，例如使用更广泛的 SIMD（AVX2、AVX512），利用连续行之间的冗余等 - 这些留给读者练习。
SSE 实现速度最快，但这是以增加复杂性、降低可维护性和降低可移植性为代价的。
OpenCV blur 函数的性能次之，如果它满足吞吐量要求，应该是首选解决方案 - 这是最简单的解决方案，简单就是好。
【讨论】：
很好的答案，但我会添加两件事。 1. 您的 SSE 版本需要 SSSE3 指令集。通常存在但尚未普遍存在，Steam 硬件调查显示“99.17%”。 2. 用于基准测试的 C++ 编译器版本和编译器选项。 @Soonts：谢谢 - 好点 - 答案已更新。以上是关于OpenCV 中的滤波函数的主要内容，如果未能解决你的问题，请参考以下文章