SSE mean filter in c++ and OpenCV

user_2152719 0
Viewed 539

I would like to modify the code for an OpenCV mean filter to use Intel intrinsics. I'm an SSE newbie and I really don't know where to start from. I checked a lot of resources on the web, but I didn't have a lot of success.

This is the program:

#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);
}

The part that I find difficult is understanding how to convert the innermost 2 loops, where the mean value is computed. Any help will be greatly appreciated.

simd sse opencv c++
1 Answers

Just for fun, I thought it might be interesting to start with a naive implementation of a 3x3 mean filter and then optimise this incrementally, ending up with a SIMD (SSE) implementation, measuring the throughput improvement at each stage.

1 - Mean_3_3_ref - reference implementation

This is just a simple scalar implementation which we'll use as a baseline for throughput and for validating further implementations:

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 - somewhat optimised scalar implementation

Exploit the redundancy in summing successive columns - we save the last two column sums so that we only need to calculate one new column sum (per channel) on each iteration:

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 - further optimised scalar implementation

As per Mean_3_3_scalar, but also remove OpenCV overheads by caching pointers to each row that we are working on:

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 - leverage OpenCV's blur function

OpenCV has a function called blur, which is based on the function boxFilter, which is just another name for a mean filter. Since OpenCV code has been quite heavily optimised over the years (using SIMD in many cases), let's see if this makes a big improvement over our scalar code:

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 implementation

This a reasonably efficient SIMD implementation. It uses the same techniques as the scalar code above in order to eliminate redundancy in processing successive pixels:

#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;
    }
}

Results

I benchmarked all of the above implementations on an 8th gen Core i9 (MacBook Pro 16,1) at 2.4 GHz, with an image size of 2337 rows x 3180 cols. The compiler was Apple clang version 12.0.5 (clang-1205.0.22.9) and the only optimisation switch was -O3. OpenCV version was 4.5.0 (via Homebrew). (Note: I verified that for Mean_3_3_blur the cv::blur function was dispatched to an AVX2 implementation.) The results:

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

Notes

  1. I have ignored the border pixels in all implementations - if required these can either be filled with pixels from the original image or using some other form of edge pixel processing.

  2. The code is not "industrial strength" - it was only written for benchmarking purposes.

  3. There are a few further possible optimisations, e.g. use wider SIMD (AVX2, AVX512), exploit the redundancy between successive rows, etc - these are left as an exercise for the reader.

  4. The SSE implementation is fastest, but this comes at the cost of increased complexity, decreased mantainability and reduced portability.

  5. The OpenCV blur function gives the second best performance, and should probably be the preferred solution if it meets throughput requirements - it's the simplest solution, and simple is good.

user_2530560
Related