Image filter for gradient detection

Question

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];
  short *gradImg = new short[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;

               gradImg[index] = sum;
               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

Answer 1

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);
    return _mm_xor_si128(_mm_add_epi16(x, m), m);
}

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);
            __m128i gx = _mm_add_epi16(_mm_add_epi16(com1, com2), _mm_sub_epi16(img5A, img6A));
            // 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 lensqL = _mm_madd_epi16(gxgyL, gxgyL);
            __m128i lensqH = _mm_madd_epi16(gxgyH, gxgyH);
            __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);

            // store gradient lengths
            size_t index = i*cols + j;
            _mm_storeu_si128((__m128i*)&gradImg[index], sum);

            // 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;

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

Answer 2

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

           gradImg[index] = sum;
           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.

Answer 3

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]();
  short *gradImg = 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;
      gradImg[i] = std::sqrt((double)(sqr));
    }
  }
}
// Register the function as a benchmark
BENCHMARK(New_Version);

This could definitely be made neater and more general (e.g. the separable convolution could be abstracted to a separate function taking the image and 2 kernels).

Note that there are also other strategies for dealing with image edges (clamp samples to the edges, mirror, wrap, zero, etc.)