C++ 和 OpenCV 中的 SSE 均值滤波器答案

【问题标题】：SSE mean filter in c++ and OpenCVC++ 和 OpenCV 中的 SSE 均值滤波器
【发布时间】：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)

标签： c++ opencv sse simd

【解决方案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：谢谢 - 好点 - 答案已更新。

1 - Mean_3_3_ref - 参考实现

2 - Mean_3_3_scalar - 一些优化的标量实现

3 - Mean_3_3_scalar_opt - 进一步优化的标量实现

4 - Mean_3_3_blur - 利用 OpenCV 的模糊功能