Быстрая реализация ковариации двух 8-битных массивов

Question

Мне нужно сравнить большое количество похожих изображений небольшого размера (до 200x200). Итак, я пытаюсь реализовать алгоритм SSIM (структурное сходство см. https://en.wikipedia.org/wiki/Structural_similarity). SSIM требует вычисления ковариации двух 8-битных серых изображений. Простейшая реализация выглядит следующим образом:

float SigmaXY(const uint8_t * x, const uint8_t * y, size_t size, float averageX, float averageY) 
{ 
    float sum = 0; 
    for(size_t i = 0; i < size; ++i) 
     sum += (x[i] - averageX) * (y[i] - averageY); 
    return sum/size; 
} 

Но он имеет низкую производительность. Поэтому я надеюсь улучшить его с помощью SIMD или CUDA (я слышал, что это можно сделать). К сожалению, у меня нет опыта для этого. Как это будет выглядеть? И куда я должен идти?

Answer 1

У меня есть другое приятное решение!

Сначала я хотел бы упомянуть некоторые математические формулы:

averageX = Sum(x[i])/size; 
averageY = Sum(y[i])/size; 

И поэтому:

Sum((x[i] - averageX)*(y[i] - averageY))/size = 

Sum(x[i]*y[i])/size - Sum(x[i]*averageY)/size - 
Sum(averageX*y[i])/size + Sum(averageX*averageY)/size = 

Sum(x[i]*y[i])/size - averageY*Sum(x[i])/size - 
averageX*Sum(y[i])/size + averageX*averageY*Sum(1)/size = 

Sum(x[i]*y[i])/size - averageY*averageX - 
averageX*averageY + averageX*averageY = 

Sum(x[i]*y[i])/size - averageY*averageX;  

Это позволяет модифицировать наш алгоритм:

float SigmaXY(const uint8_t * x, const uint8_t * y, size_t size, float averageX, float averageY) 
{ 
    uint32_t sum = 0; // If images will have size greater then 256x256 than you have to use uint64_t. 
    for(size_t i = 0; i < size; ++i) 
     sum += x[i]*y[i]; 
    return sum/size - averageY*averageX; 
} 

И только после этого мы может использовать SIMD (я использовал SSE2):

#include <emmintrin.h> 

inline __m128i SigmaXY(__m128i x, __m128i y) 
{ 
    __m128i lo = _mm_madd_epi16(_mm_unpacklo_epi8(x, _mm_setzero_si128()), _mm_unpacklo_epi8(y, _mm_setzero_si128())); 
    __m128i hi = _mm_madd_epi16(_mm_unpackhi_epi8(y, _mm_setzero_si128()), _mm_unpackhi_epi8(y, _mm_setzero_si128())); 
    return _mm_add_epi32(lo, hi); 
} 

float SigmaXY(const uint8_t * x, const uint8_t * y, size_t size, float averageX, float averageY) 
{ 
    uint32_t sum = 0; 
    size_t i = 0, alignedSize = size/16*16; 
    if(size >= 16) 
    { 
     __m128i sums = _mm_setzero_si128(); 
     for(; i < alignedSize; i += 16) 
     { 
      __m128i _x = _mm_loadu_si128((__m128i*)(x + i)); 
      __m128i _y = _mm_loadu_si128((__m128i*)(y + i)); 
      sums = _mm_add_epi32(sums, SigmaXY(_x, _y)); 
     } 
     uint32_t _sums[4]; 
     _mm_storeu_si128(_sums, sums); 
     sum = _sums[0] + _sums[1] + _sums[2] + _sums[3]; 
    } 
    for(; i < size; ++i) 
     sum += x[i]*y[i]; 
    return sum/size - averageY*averageX; 
} 

Answer 2

Существует реализация SIMD алгоритма (я использовал SSE4.1):

#include <smmintrin.h> 

template <int shift> inline __m128 SigmaXY(const __m128i & x, const __m128i & y, __m128 & averageX, __m128 & averageY) 
{ 
    __m128 _x = _mm_cvtepi32_ps(_mm_cvtepu8_epi32(_mm_srli_si128(x, shift))); 
    __m128 _y = _mm_cvtepi32_ps(_mm_cvtepu8_epi32(_mm_srli_si128(y, shift))); 
    return _mm_mul_ps(_mm_sub_ps(_x, averageX), _mm_sub_ps(_y, averageY)) 
} 

float SigmaXY(const uint8_t * x, const uint8_t * y, size_t size, float averageX, float averageY) 
{ 
    float sum = 0; 
    size_t i = 0, alignedSize = size/16*16; 
    if(size >= 16) 
    { 
     __m128 sums = _mm_setzero_ps(); 
     __m128 avgX = _mm_set1_ps(averageX); 
     __m128 avgY = _mm_set1_ps(averageY); 
     for(; i < alignedSize; i += 16) 
     { 
      __m128i _x = _mm_loadu_si128((__m128i*)(x + i)); 
      __m128i _y = _mm_loadu_si128((__m128i*)(y + i)); 
      sums = _mm_add_ps(sums, SigmaXY<0>(_x, _y, avgX, avgY); 
      sums = _mm_add_ps(sums, SigmaXY<4>(_x, _y, avgX, avgY); 
      sums = _mm_add_ps(sums, SigmaXY<8>(_x, _y, avgX, avgY); 
      sums = _mm_add_ps(sums, SigmaXY<12>(_x, _y, avgX, avgY); 
     } 
     float _sums[4]; 
     _mm_storeu_ps(_sums, sums); 
     sum = _sums[0] + _sums[1] + _sums[2] + _sums[3]; 
    } 
    for(; i < size; ++i) 
     sum += (x[i] - averageX) * (y[i] - averageY); 
    return sum/size; 
} 

Я надеюсь, что это будет полезно для вас.