December 3, 2013 at 17:26Optimization of image processing in C ++ using SIMD. Median filter
- Tutorial
Introduction
Earlier in the introductory articleI raised a list of problems that a developer will have to face if he wants to optimize image processing optimization using SIMD instructions. Now the time has come to show by a concrete example how the above problems can be solved. I thought for a long time which algorithm to choose for the first example, and decided to focus on median filtering. Median filtering is an effective way to suppress noise that inevitably appears on digital cameras in low-light scenes. This algorithm is quite resource-intensive - for example, when processing a gray image with a 3x3 median filter, about 50 operations per image point are required. But at the same time, he only operates with 8-bit numbers and he needs relatively little input to work.
For reference, I recall the essence of the algorithm:
- For each point of the original image, a certain neighborhood is taken (in our case, 3x3).
- Points in a given neighborhood are sorted by increasing brightness.
- The midpoint (5th for a 3x3 filter) of the sorted neighborhood is recorded in the final image.
For the first example, I decided to simplify the task and consider the case of a gray image (8 bits per point) without taking into account boundary effects - for this, the central part of a slightly larger image was taken for filtering. Of course, in real problems, the boundary points of the image should be calculated using a special algorithm, and sometimes the implementation of these special algorithms can take up most of the code, but we will dwell on these issues another time. In addition, for similar reasons, it is assumed that the image width is a multiple of 32.
Scalar versions of the algorithm
1st scalar version
In this initial implementation, the problem was solved head on: the neighborhood for each point of the original image was copied into an auxiliary array, which was then sorted using the std :: sort () function.
inline void Load(const uint8_t * src, int a[3])
{
a[0] = src[-1];
a[1] = src[0];
a[2] = src[1];
}
inline void Load(const uint8_t * src, size_t stride, int a[9])
{
Load(src - stride, a + 0);
Load(src, a + 3);
Load(src + stride, a + 6);
}
inline void Sort(int a[9])
{
std::sort(a, a + 9);
}
void MedianFilter(const uint8_t * src, size_t srcStride, size_t width, size_t height, uint8_t * dst, size_t dstStride)
{
int a[9];
for(size_t y = 0; y < height; ++y)
{
for(size_t x = 0; x < width; ++x)
{
Load(src + x, srcStride, a);
Sort(a);
dst[x] = (uint8_t)a[4];
}
src += srcStride;
dst += dstStride;
}
}
2nd scalar version
The previous method, however, immediately shows a bottleneck - this is the standard sorting function. It is designed to sort an array of arbitrary size, it performs many internal checks, so it most likely will not be the optimal solution for this problem, where the array is small and has a fixed size. In addition, it is not at all necessary for us to completely sort the array, it is enough that the value we want is in the 4th element. Therefore, we can apply the “network” sorting method (also known as the bitonic sorting method):
inline void Sort(int & a, int & b) // сортирует пару чисел
{
if(a > b)
{
int t = a;
a = b;
b = t;
}
}
inline void Sort(int a[9]) //частично сортирует весь массив
{
Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]);
Sort(a[0], a[1]); Sort(a[3], a[4]); Sort(a[6], a[7]);
Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]);
Sort(a[0], a[3]); Sort(a[5], a[8]); Sort(a[4], a[7]);
Sort(a[3], a[6]); Sort(a[1], a[4]); Sort(a[2], a[5]);
Sort(a[4], a[7]); Sort(a[4], a[2]); Sort(a[6], a[4]);
Sort(a[4], a[2]);
}
3rd scalar version
Although the first method turned out to be much faster than the zero method (see the test result below), it still has one bottleneck - there are a lot of conditional transitions in this method. As you know, modern processors do not like them very much, because they have a pipelined architecture and the possibility of the extraordinary execution of several instructions in one clock cycle. For both the first and the second, conditional transitions are extremely undesirable, since the processor must stop and wait for the results of calculations in order to find out which branch of the program it should go next. Fortunately, sorting two 8-bit values can be implemented without conditional branching either:
inline void Sort(int & a, int & b)
{
int d = a - b;
int m = ~(d >> 8);
b += d&m;
a -= d&m;
}
The result of these preliminary optimizations of the scalar code is the following version of the median filter:
inline void Load(const uint8_t * src, int a[3])
{
a[0] = src[-1];
a[1] = src[0];
a[2] = src[1];
}
inline void Load(const uint8_t * src, size_t stride, int a[9])
{
Load(src - stride, a + 0);
Load(src, a + 3);
Load(src + stride, a + 6);
}
inline void Sort(int & a, int & b)
{
int d = a - b;
int m = ~(d >> 8);
b += d&m;
a -= d&m;
}
inline void Sort(int a[9]) //частично сортирует весь массив
{
Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]);
Sort(a[0], a[1]); Sort(a[3], a[4]); Sort(a[6], a[7]);
Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]);
Sort(a[0], a[3]); Sort(a[5], a[8]); Sort(a[4], a[7]);
Sort(a[3], a[6]); Sort(a[1], a[4]); Sort(a[2], a[5]);
Sort(a[4], a[7]); Sort(a[4], a[2]); Sort(a[6], a[4]);
Sort(a[4], a[2]);
}
void MedianFilter(const uint8_t * src, size_t srcStride, size_t width, size_t height, uint8_t * dst, size_t dstStride)
{
int a[9];
for(size_t y = 0; y < height; ++y)
{
for(size_t x = 0; x < width; ++x)
{
Load(src + x, srcStride, a);
Sort(a);
dst[x] = (uint8_t)a[4];
}
src += srcStride;
dst += dstStride;
}
}
This version is already significantly (somewhere in 5-6 times) faster than our initial version of the algorithm in speed. And it is based on it that we will carry out SIMD optimization of the algorithm, as well as compare the speed of the scalar and vector versions of the algorithm.
SIMD version of the algorithm
To illustrate the optimization of the median filtering algorithm using SIMD, I will use two sets of instructions SSE2 and AVX2, missing the MMX extension, which is currently outdated and has more historical interest. Fortunately, in order to use SIMD instructions, it is not necessary to use assembler . Most modern C ++ compilers have support for intrinsics (built-in functions with which you can use all kinds of processor extensions). Programming using intrinsics is practically no different from programming on pure C. Intrinsic functions are converted by the compiler directly into processor instructions, although working directly with processor registers remains hidden from the programmer.
SSE2 version
SSE2 integer commands are defined in the header file <emmintrin.h>. The basic type is __m128i - a 128 bit vector, which, depending on the context, can be interpreted as a set of 2 64 bit, 4 32 bit, 8 x 16 bit, 16 x 8 bit signed or unsigned numbers. As you can see, they support not only vector arithmetic operations, but also vector logical operations, as well as vector operations of loading and unloading data. The following are optimizations for the median filter using SSE2 instructions. The code, it seems to me, is pretty visual.
#include < emmintrin.h >
inline void Load(const uint8_t * src, __m128i a[3])
{
a[0] = _mm_loadu_si128((__m128i*)(src - 1)); //загрузка 128 битного вектора по невыровненному по 16 битной границе адресу
a[1] = _mm_loadu_si128((__m128i*)(src));
a[2] = _mm_loadu_si128((__m128i*)(src + 1));
}
inline void Load(const uint8_t * src, size_t stride, __m128i a[9])
{
Load(src - stride, a + 0);
Load(src, a + 3);
Load(src + stride, a + 6);
}
inline void Sort(__m128i& a, __m128i& b)
{
__m128i t = a;
a = _mm_min_epu8(t, b); //нахождение минимума 2-х 8 битных беззнаковых чисел для каждого из 16 значений вектора
b = _mm_max_epu8(t, b); //нахождение максимума 2-х 8 битных беззнаковых чисел для каждого из 16 значений вектора
}
inline void Sort(__m128i a[9]) //частично сортирует весь массив
{
Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]);
Sort(a[0], a[1]); Sort(a[3], a[4]); Sort(a[6], a[7]);
Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]);
Sort(a[0], a[3]); Sort(a[5], a[8]); Sort(a[4], a[7]);
Sort(a[3], a[6]); Sort(a[1], a[4]); Sort(a[2], a[5]);
Sort(a[4], a[7]); Sort(a[4], a[2]); Sort(a[6], a[4]);
Sort(a[4], a[2]);
}
void MedianFilter(const uint8_t * src, size_t srcStride, size_t width, size_t height, uint8_t * dst, size_t dstStride)
{
__m128i a[9];
for(size_t y = 0; y < height; ++y)
{
for(size_t x = 0; x < width; x += sizeof(__m128i))
{
Load(src + x, srcStride, a);
Sort(a);
_mm_storeu_si128((__m128i*)(dst + x), a[4]); //сохранение 128 битного вектора по невыровненному по 16 битной границе адресу
}
src += srcStride;
dst += dstStride;
}
}
AVX2 version
The integer AVX2 commands are defined in the header file <immintrin.h>. The basic type is __m256i - a 256 bit vector, which, depending on the context, can be interpreted as a set of 4 64 bit, 8 x 32 bit, 16 x 16 bit, 32 x 8 bit signed or unsigned numbers. Although the AVX2 instruction set largely repeats the SSE2 instruction set (given the doubled width of the vector), it also contains quite a few instructions that have no analogues in SSE2. Below is the optimization of the median filter using the AVX2 instructions.
#include < immintrin.h >
inline void Load(const uint8_t * src, __m256i a[3])
{
a[0] = _mm256_loadu_si256((__m128i*)(src - 1)); //загрузка 256 битного вектора по невыровненному по 32 битной границе адресу
a[1] = _mm256_loadu_si256((__m128i*)(src));
a[2] = _mm256_loadu_si256((__m128i*)(src + 1));
}
inline void Load(const uint8_t * src, size_t stride, __m256i a[9])
{
Load(src - stride, a + 0);
Load(src, a + 3);
Load(src + stride, a + 6);
}
inline void Sort(__m256i& a, __m256i& b)
{
__m256i t = a;
a = _mm256_min_epu8(t, b); //нахождение минимума 2-х 8 битных беззнаковых чисел для каждого из 32 значений вектора
b = _mm256_max_epu8(t, b); //нахождение максимума 2-х 8 битных беззнаковых чисел для каждого из 32 значений вектора
}
inline void Sort(__m256i a[9]) //частично сортирует весь массив
{
Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]);
Sort(a[0], a[1]); Sort(a[3], a[4]); Sort(a[6], a[7]);
Sort(a[1], a[2]); Sort(a[4], a[5]); Sort(a[7], a[8]);
Sort(a[0], a[3]); Sort(a[5], a[8]); Sort(a[4], a[7]);
Sort(a[3], a[6]); Sort(a[1], a[4]); Sort(a[2], a[5]);
Sort(a[4], a[7]); Sort(a[4], a[2]); Sort(a[6], a[4]);
Sort(a[4], a[2]);
}
void MedianFilter(const uint8_t * src, size_t srcStride, size_t width, size_t height, uint8_t * dst, size_t dstStride)
{
__m256i a[9];
for(size_t y = 0; y < height; ++y)
{
for(size_t x = 0; x < width; x += sizeof(__m256i))
{
Load(src + x, srcStride, a);
Sort(a);
_mm256_storeu_si256((__m256i*)(dst + x), a[4]); //сохранение 256 битного вектора по невыровненному по 32 битной границе адресу
}
src += srcStride;
dst += dstStride;
}
}
Note: To get the maximum acceleration from optimization, code that contains AVX2 instructions must be compiled with the following compiler options: (/ arch: AVX for Visual Studio 2012, -march = core-avx2 for GCC). In addition, it is very desirable that the code does not contain alternating AVX2 and SSE2 instructions, since switching the processor operating mode in this case will adversely affect overall performance. Based on the foregoing, it is advisable to place the AVX2 and SSE2 versions of the algorithms in different ".cpp" files.
Testing
Testing was performed on 2 MB gray images (1920 x 1080). To improve the accuracy of measuring time, tests were run several times. The execution time was obtained by dividing the total test execution time by the number of runs. The number of runs was chosen in such a way that the total execution time was not less than 1 second, which should provide two signs of accuracy when measuring the execution time of the algorithms. The algorithms were compiled with maximum optimization for 64 bit Windows 7 on Microsoft Visual Studio 2012 and run on the iCore-7 4770 processor (3.4 GHz).
| Execution time, ms | Relative acceleration | ||||
| Best scalar | SSE2 | AVX2 | Best scalar / SSE2 | Best scalar / AVX2 | SSE2 / AVX2 |
| 24.814 | 0.565 | 0.424 | 43.920 | 58.566 | 1.333 |
Conclusion
As can be seen from the test results, the acceleration of the median filtering algorithm from the use of SIMD instructions on modern processors can reach 40 to 60 times. It seems to me that this is a sufficient amount to bother with optimization (unless of course execution speed is important in your program). In this publication, I tried to simplify the code used as much as possible.
So, issues related to data alignment, processing of boundary points, and much more remained beyond the scope.
I will try to disclose these questions in the following articles. If the reader is interested in what the battle code that implements this functionality will look like, he will be able to find it here .
