# Image filter for gradient detection

This filter is basically gradient detection (the variable gradImg) like an edge detection (Sobel or Prewitt filters) and also the direction of that gradient whether vertical or horizontal (dirImg).

I'm trying to optimize it to run faster and I found that this part of the code slow. Would you have any suggestions please? Here's the code you can copy to http://quick-bench.com/.

#include <math.h>
static void Old_Version(benchmark::State& state) {
// Code inside this loop is measured repeatedly
int rows = 1000, cols = 1000;
unsigned char *smoothImg = new unsigned char[cols*rows];
unsigned char* dirImg = new unsigned char[cols*rows];
for (auto _ : state) {
for (int i = 1; i<rows - 1; i++) {
for (int j = 1; j<cols - 1; j++) {
int com1 = smoothImg[(i + 1)*cols + j + 1] - smoothImg[(i - 1)*cols + j - 1];
int com2 = smoothImg[(i - 1)*cols + j + 1] - smoothImg[(i + 1)*cols + j - 1];

int gx=abs(com1+com2+(smoothImg[i*cols + j + 1] - smoothImg[i*cols + j - 1]));
int gy=abs(com1-com2+(smoothImg[(i + 1)*cols + j] - smoothImg[(i - 1)*cols + j]));

int sum = (int)sqrt((double)gx*gx + gy*gy);

int index = i*cols + j;

if (sum >= 20) {
if (gx >= gy) dirImg[index] = 1;//1 vertical
else      dirImg[index] = 2;//2 Horizontal
}
}
}
}
}
// Register the function as a benchmark
BENCHMARK(Old_Version);


My compilation command, using GCC (Ubuntu 8.1.0-5ubuntu1~16.04) 8.1.0:

g++ -std=c++1z -fomit-frame-pointer -O4 -ffast-math -mmmx -msse -msse2 -msse3 -DNDEBUG -Wall improve_code.cpp -o improve_code -fopenmp

• Could you add some details about what this filter does? Feb 14, 2019 at 14:16
• Looks like smoothImg contents are used uninitialized, which makes for Undefined behaviour. Is this the real code from your project? Feb 14, 2019 at 14:19
• @user673679 I added some information Feb 14, 2019 at 14:21
• @TobySpeight you're right, but smoothImg in my code is an image .. you can imagine a random values Feb 14, 2019 at 14:22
• Are SSE/AVX intrinsics within the scope of this question? Feb 14, 2019 at 14:48

For SIMD optimization, I tried doing it with SSE2, as the compilation flags imply SSSE3 is disabled, and SSE3 is not useful here.

It's mostly a transliteration of the scalar code, but there are a few points of interest:

• Except the squaring of gx and gy which is done in 32bit, and taking the square root which is done in floating point, most arithmetic is done in 16bit. That shouldn't change the results as everything fits well within the margins.
• 16bit is more than 8bit, so basically there is a choice between doing some half-wide loads and stores, or "doubling" the code. Often "doubling" the code is a little faster because it gets more arithmetic done with the same number of loads and stores, but that didn't pan out this time. So I used the less-common _mm_loadl_epi64 and the corresponding store. Note that these notionally take a pointer to __m128i, but they actually just load and store a qword (with no alignment restriction).
• _mm_madd_epi16 is used for the gx*gx + gy*gy expression, which requires interleaving gx and gy, but is still more convenient than the alternatives.
• There is no good way to conditionally not write a byte somewhere in the middle of a qword store, so for gradients less than 20 I write a zero, which is definitely different than the scalar code (but maybe safer anyway)
• SSE2 has no single instruction to take an absolute value, but SSSE3 does. I didn't use it, but I factored it out so it can easily be updated if SSSE3 support can be assumed someday.
• I paid no particular attention to alignment. The qword loads and stores don't care that much about alignment inherently, but on Core2 it is nevertheless Quite Bad to cross a cache line that way. If Core2 is a serious target, perhaps that aspect can be improved.
• No testing was done other than for performance.

The results on quick-bench are OK but not amazing, about an 4x improvement (the final code there is labeled Gradient2). But maybe that's all we can hope for, given that there are square roots to be done. Of course you can still add threading on top of this which should (mostly) stack multiplicatively.

#include <x86intrin.h>

__m128i abs_epi16(__m128i x)
{
__m128i m = _mm_srai_epi16(x, 15);
}

void Gradient(unsigned char *smoothImg, short *gradImg, unsigned char* dirImg, size_t rows, size_t cols) {
for (size_t i = 1; i + 1 < rows; i++) {
size_t j = 1;
// do blocks of 8 pixels at the time until near the edge
for (; j + 8 < cols; j += 8) {
__m128i zero = _mm_setzero_si128();
// com1
__m128i img1 = _mm_loadl_epi64((__m128i*)&smoothImg[(i + 1)*cols + j + 1]);
__m128i img2 = _mm_loadl_epi64((__m128i*)&smoothImg[(i - 1)*cols + j - 1]);
__m128i img1A = _mm_unpacklo_epi8(img1, zero);
__m128i img2A = _mm_unpacklo_epi8(img2, zero);
__m128i com1 = _mm_sub_epi16(img1A, img2A);
// com2
__m128i img3 = _mm_loadl_epi64((__m128i*)&smoothImg[(i - 1)*cols + j + 1]);
__m128i img4 = _mm_loadl_epi64((__m128i*)&smoothImg[(i + 1)*cols + j - 1]);
__m128i img3A = _mm_unpacklo_epi8(img3, zero);
__m128i img4A = _mm_unpacklo_epi8(img4, zero);
__m128i com2 = _mm_sub_epi16(img3A, img4A);
// gx
__m128i img5 = _mm_loadl_epi64((__m128i*)&smoothImg[i*cols + j + 1]);
__m128i img6 = _mm_loadl_epi64((__m128i*)&smoothImg[i*cols + j - 1]);
__m128i img5A = _mm_unpacklo_epi8(img5, zero);
__m128i img6A = _mm_unpacklo_epi8(img6, zero);
// gy
__m128i img7 = _mm_loadl_epi64((__m128i*)&smoothImg[(i + 1)*cols + j]);
__m128i img8 = _mm_loadl_epi64((__m128i*)&smoothImg[(i - 1)*cols + j]);
__m128i img7A = _mm_unpacklo_epi8(img7, zero);
__m128i img8A = _mm_unpacklo_epi8(img8, zero);
__m128i gy = _mm_add_epi16(_mm_sub_epi16(com1, com2), _mm_sub_epi16(img7A, img8A));
// sum
// gx and gy are interleaved, multiplied by themselves and
// horizontally added in pairs, creating gx*gx+gy*gy as a dword
// 32bits is required here to avoid overflow, but also convenient for the next step
__m128i gxgyL = _mm_unpacklo_epi16(gx, gy);
__m128i gxgyH = _mm_unpackhi_epi16(gx, gy);
__m128i lenL = _mm_cvttps_epi32(_mm_sqrt_ps(_mm_cvtepi32_ps(lensqL)));
__m128i lenH = _mm_cvttps_epi32(_mm_sqrt_ps(_mm_cvtepi32_ps(lensqH)));
__m128i sum = _mm_packs_epi32(lenL, lenH);

size_t index = i*cols + j;

// classify H/V/low
__m128i thresholdLow = _mm_set1_epi16(19);
__m128i markerV = _mm_set1_epi8(1);
__m128i isSignificant = _mm_cmpgt_epi16(sum, thresholdLow);
__m128i isHorizontal = _mm_cmplt_epi16(abs_epi16(gx), abs_epi16(gy));
// if not horizontal, then 1 - 0 = 1
// if horizontal,  then 1 - (-1) = 2
// if not significant, then make it zero
__m128i classifier = _mm_and_si128(_mm_sub_epi16(markerV, isHorizontal), isSignificant);
_mm_storel_epi64((__m128i*)&dirImg[index], _mm_packs_epi16(classifier, classifier));
}
for (; j + 1 < cols; j++) {
int com1 = smoothImg[(i + 1)*cols + j + 1] - smoothImg[(i - 1)*cols + j - 1];
int com2 = smoothImg[(i - 1)*cols + j + 1] - smoothImg[(i + 1)*cols + j - 1];

int gx=abs(com1+com2+(smoothImg[i*cols + j + 1] - smoothImg[i*cols + j - 1]));
int gy=abs(com1-com2+(smoothImg[(i + 1)*cols + j] - smoothImg[(i - 1)*cols + j]));

int sum = (int)sqrt((double)gx*gx + gy*gy);

size_t index = i*cols + j;

if (sum >= 20) {
if (gx >= gy) dirImg[index] = 1; //1 vertical
else      dirImg[index] = 2; //2 Horizontal
}
}
}
}

• Thank you for your time.. that's what I'm looking for.. It's faster I just need to test it on real image. I need to learn more about SIMD optimization, it is very useful Feb 15, 2019 at 11:17

Prefer C++ headers (<cmath> rather than <math.h>); we'll then use std::abs() and std::sqrt() - or better, std::hypot().

The std::abs() calls seem premature - we can defer them to within the (sum >= 20) condition.

Prefer to give names to the values used in dirImg. For example:

enum orientation : unsigned char
{
None,
Vertical,
Horizontal,
};


The indexing in the loop might be easier to read if we just offset from index; it might also be slightly more efficient if the compiler can't reason about the (cols ± 1) * i) for us:

  for (std::size_t i = 1; i<rows - 1; i++) {
for (std::size_t j = 1; j<cols - 1; j++) {
auto const index = i*cols + j;

// leading and trailing diagonal differences are common
int com1 = smoothImg[index + cols + 1] - smoothImg[index - cols - 1];
int com2 = smoothImg[index - cols + 1] - smoothImg[index + cols - 1];

int gx = com1 + com2 + smoothImg[index + 1] - smoothImg[index - 1];
int gy = com1 - com2 + smoothImg[index + cols] - smoothImg[index - cols];

auto sum = static_cast<short>(std::hypot(gx, gy));

if (sum >= 20) {
dirImg[index] = std::abs(gx) >= std::abs(gy) ? Vertical : Horizontal;
}
}
}


What's special about the threshold value 20? Perhaps that should be a parameter to the function.

Remember to delete[] what you new[] - or better, use standard containers or smart pointers so that we don't need to remember, and so the storage is reclaimed even on a non-local exit (e.g. due to a later std::bad_alloc).

Consider spreading the work across processor cores, e.g. by applying #pragma omp parallel for to the i loop.

Be careful if you use the outputs - we've left a lot of the values uninitialised, especially in dirImg.

• Is hypot really supposed to be better? A quick test suggests it's substantially worse Feb 14, 2019 at 17:18
• Yes, std::hypot() is required to be accurate to within 1 ULP, and there are tighter constrains on returning subnormal values. It's also easier to read. Feb 14, 2019 at 17:22
• But it makes the whole thing 2.5 as slow. That doesn't work on a question that asks how to make it faster. Besides, the result is truncated and the input cannot be subnormal (apart from zero but do we really consider that a subnormal?) Feb 14, 2019 at 17:41
• @harold - the question wasn't tagged performance until after I published this answer. Feb 14, 2019 at 17:42
• It's more correct, but you always have the choice of getting incorrect answers more quickly... Feb 15, 2019 at 10:29

Both Prewitt and Sobel filters are "separable".

To calculate the gradients (gx and gy), the original code effectively iterates over and sums a 2d kernel (i.e. convolution). For the Prewitt and Sobel filters, each gradient kernel can be split into two 1D kernels, and we can do the calculation of each gradient in two passes. For a k sized filter, this reduces the work needed at each pixel from k * k to 2 * k.

This does require some extra memory, but the following is around 4.8x faster.

static void New_Version(benchmark::State& state) {

const int rows = 1000;
const int cols = 1000;
const int threshold = 20;

auto to_index = [&] (std::size_t x, std::size_t y)
{
return y * cols + x;
};

unsigned char *smoothImg = new unsigned char[cols*rows];
short *scratch = new short[cols*rows](); // note: could just use gradImg instead.
short *gx = new short[cols*rows]();
short *gy = new short[cols*rows]();
unsigned char* dirImg = new unsigned char[cols*rows];

const short gxx[3] = { +1,  0, -1 };
const short gxy[3] = { +1, +1, +1 };

const short gyx[3] = { +1, +1, +1 };
const short gyy[3] = { +1,  0, -1 };

for (auto _ : state) {

// x gradient: convolve in the x direction first
for (std::size_t y = 0; y != rows; ++y) // note: need to calculate for the edge pixels too
for (std::size_t x = 1; x != cols - 1; ++x)
for (std::size_t o = 0; o != 3; ++o)
scratch[to_index(x, y)] += smoothImg[to_index(x + o - 1, y)] * gxx[o];

// x gradient: use the results of the first pass and convolve in the y direction
for (std::size_t y = 1; y != rows - 1; ++y)
for (std::size_t x = 0; x != cols; ++x)
for (std::size_t o = 0; o != 3; ++o)
gx[to_index(x, y)] += scratch[to_index(x, y + o - 1)] * gxy[o];

// y gradient: convolve in the x direction first
for (std::size_t y = 0; y != rows; ++y)
for (std::size_t x = 1; x != cols - 1; ++x)
{
scratch[to_index(x, y)] = 0;
for (std::size_t o = 0; o != 3; ++o)
scratch[to_index(x, y)] += smoothImg[to_index(x + o - 1, y)] * gyx[o];
}

// y gradient: use the results of the first pass and convolve in the y direction
for (std::size_t y = 1; y != rows - 1; ++y)
for (std::size_t x = 0; x != cols; ++x)
for (std::size_t o = 0; o != 3; ++o)
gy[to_index(x, y)] += scratch[to_index(x, y + o - 1)] * gyy[o];

// calculate magnitude and direction:
const unsigned char v = 1;
const unsigned char h = 2;
const int t2 = threshold * threshold;
for (std::size_t i = 0; i != cols * rows; ++i)
{
auto const x2 = gx[i] * gx[i];
auto const y2 = gy[i] * gy[i];
auto const sqr = x2 + y2;
dirImg[i] = (sqr >= t2) ? (x2 >= y2) ? v : h : 0;