c++ - C++ 和 OpenCV 中的 SSE 均值滤波器

Question

我想修改 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 个循环，其中计算平均值。任何帮助将不胜感激。

Answer 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 - 一些优化的标量实现
利用对连续列求和的冗余——我们保存最后两列的总和，以便我们只需要在每次迭代中计算一个新的列总和(每个 channel ):

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 函数已分派(dispatch)到 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 function 提供次佳的性能，如果它满足吞吐量要求，它可能应该是首选解决方案 - 这是最简单的解决方案，简单就是好。