Быстрый расчет моментов изображения

Question

У меня есть маска (8-битное серое изображение), и мне нужно вычислить центр области с заданным индексом маски. Для этого мне нужно рассчитать моменты первого порядка по осям X и Y для этой маски. В настоящее время я использую следующий код:

void GetCenter(const uint8_t * mask, size_t stride, size_t width, size_t height, 
       uint8_t index, double * centerX, double * centerY) 
{ 
    uint64_t sum = 0, sumX = 0, sumY = 0; 
    for(size_t y = 0; y < height; ++y) 
    { 
     for(size_t x = 0; x < width; ++x) 
     { 
      if(mask[x] == index) 
      { 
       sum++; 
       sumX += x; 
       sumY += y; 
      }    
     } 
     mask += stride; 
    } 
    *centerX = sum ? (double)sumX/sum : 0; 
    *centerY = sum ? (double)sumY/sum : 0; 
} 

И у меня вопрос: Есть ли способ, чтобы улучшить производительность этого алгоритма?

Answer 1

Существует способ значительно (более десяти раз) улучшить производительность этого алгоритма. Для этого вам необходимо использовать SIMD-инструкции для процессора, такие как (SSE2, AVX2, Altivec, NEON и т. Д.). я написал пример с использованием инструкций SSE2 (AVX2 код будет похож на него):

const __m128i K_0 = _mm_setzero_si128(); 
const __m128i K8_1 = _mm_set1_epi8(1); 
const __m128i K16_1 = _mm_set1_epi16(1); 
const __m128i K16_8 = _mm_set1_epi16(8); 
const __m128i K16_I = _mm_setr_epi16(0, 1, 2, 3, 4, 5, 6, 7); 

inline void AddMoments(const __m128i & mask, const __m128i & x, const __m128i & y, 
    __m128i & sumX, __m128i & sumY) 
{ 
    sumX = _mm_add_epi32(sumX, _mm_madd_epi16(_mm_and_si128(mask, x), K16_1)); 
    sumY = _mm_add_epi32(sumY, _mm_madd_epi16(_mm_and_si128(mask, y), K16_1)); 
} 

inline int ExtractSum(__m128i a) 
{ 
    return _mm_cvtsi128_si32(a) + _mm_cvtsi128_si32(_mm_srli_si128(a, 4)) + 
     _mm_cvtsi128_si32(_mm_srli_si128(a, 8)) + _mm_cvtsi128_si32(_mm_srli_si128(a, 12)); 
} 

void GetCenter(const uint8_t * mask, size_t stride, size_t width, size_t height, 
    uint8_t index, double * centerX, double * centerY) 
{ 
    size_t alignedWidth = width & ~(sizeof(__m128i) - 1); 
    const __m128i _index = _mm_set1_epi8(index); 

    uint64_t sum = 0, sumX = 0, sumY = 0; 
    for(size_t y = 0; y < height; ++y) 
    { 
     size_t x = 0; 

     __m128i _x = K16_I; 
     __m128i _y = _mm_set1_epi16((short)y); 
     __m128i _sum = K_0; 
     __m128i _sumX = K_0; 
     __m128i _sumY = K_0; 
     for(; x < alignedWidth; x += sizeof(__m128i)) 
     { 
      __m128i _mask = _mm_and_si128(_mm_cmpeq_epi8(_mm_loadu_si128((__m128i*)(mask + x)), _index), K8_1); 
      _sum = _mm_add_epi64(_sum, _mm_sad_epu8(_mask, K_0)); 
      AddMoments(_mm_cmpeq_epi16(_mm_unpacklo_epi8(_mask, K_0), K16_1), _x, _y, _sumX, _sumY); 
      _x = _mm_add_epi16(_x, K16_8); 
      AddMoments(_mm_cmpeq_epi16(_mm_unpackhi_epi8(_mask, K_0), K16_1), _x, _y, _sumX, _sumY); 
      _x = _mm_add_epi16(_x, K16_8); 
     } 
     sum += ExtractSum(_sum); 
     sumX += ExtractSum(_sumX); 
     sumY += ExtractSum(_sumY); 

     for(; x < width; ++x) 
     { 
      if(mask[x] == index) 
      { 
       sum++; 
       sumX += x; 
       sumY += y; 
      }    
     } 
     mask += stride; 
    } 
    *centerX = sum ? (double)sumX/sum : 0; 
    *centerY = sum ? (double)sumY/sum : 0; 
} 

P.S. Существует более простой и кросс-платформенный способ повышения производительности с использованием внешней библиотеки (http://simd.sourceforge.net/):

void GetCenter(const uint8_t * mask, size_t stride, size_t width, size_t height, 
    uint8_t index, double * centerX, double * centerY) 
{ 
    uint64_t sum, sumX, sumY, sumXX, sumXY, sumYY; 
    ::SimdGetMoments(mask, stride, width, height, index, 
     &sum, &sumX, &sumY, &sumXX, &sumXY, &sumYY); 
    *centerX = sum ? (double)sumX/sum : 0; 
    *centerY = sum ? (double)sumY/sum : 0; 
} 

Реализация с использованием _mm_movemask_epi8 и 8-битные таблиц поиска:

uint8_t g_sum[1 << 8], g_sumX[1 << 8]; 
bool Init() 
{ 
    for(int i = 0, n = 1 << 8; i < n; ++i) 
    { 
     g_sum[i] = 0; 
     g_sumX[i] = 0; 
     for(int j = 0; j < 8; ++j) 
     { 
      g_sum[i] += (i >> j) & 1; 
      g_sumX[i] += ((i >> j) & 1)*j; 
     } 
    } 
    return true; 
} 
bool g_inited = Init(); 

inline void AddMoments(uint8_t mask, size_t x, size_t y, 
         uint64_t & sum, uint64_t & sumX, uint64_t & sumY) 
{ 
    int value = g_sum[mask]; 
    sum += value; 
    sumX += x * value + g_sumX[mask]; 
    sumY += y * value; 
} 

void GetCenter(const uint8_t * mask, size_t stride, size_t width, size_t height, 
    uint8_t index, double * centerX, double * centerY) 
{ 
    size_t alignedWidth = width & ~(sizeof(__m128i) - 1); 
    const __m128i _index = _mm_set1_epi8(index); 

    union PackedValue 
    { 
     uint8_t u8[4]; 
     uint16_t u16[2]; 
     uint32_t u32; 
    } _mask; 

    uint64_t sum = 0, sumX = 0, sumY = 0; 
    for(size_t y = 0; y < height; ++y) 
    { 
     size_t x = 0; 

     for(; x < alignedWidth; x += sizeof(__m128i)) 
     { 
      _mask.u32 = _mm_movemask_epi8(_mm_cmpeq_epi8(
         _mm_loadu_si128((__m128i*)(mask + x)), _index)); 
      AddMoments(_mask.u8[0], x, y, sum, sumX, sumY); 
      AddMoments(_mask.u8[1], x + 8, y, sum, sumX, sumY); 
     } 

     for(; x < width; ++x) 
     { 
      if(mask[x] == index) 
      { 
       sum++; 
       sumX += x; 
       sumY += y; 
      }    
     } 
     mask += stride; 
    } 
    *centerX = sum ? (double)sumX/sum : 0; 
    *centerY = sum ? (double)sumY/sum : 0; 
} 

Реализация с использованием _mm_movemask_epi8 и 16-разрядных таблиц поиска:

uint16_t g_sum[1 << 16], g_sumX[1 << 16]; 
bool Init() 
{ 
    for(int i = 0, n = 1 << 16; i < n; ++i) 
    { 
     g_sum[i] = 0; 
     g_sumX[i] = 0; 
     for(int j = 0; j < 16; ++j) 
     { 
      g_sum[i] += (i >> j) & 1; 
      g_sumX[i] += ((i >> j) & 1)*j; 
     } 
    } 
    return true; 
} 
bool g_inited = Init(); 

inline void AddMoments(uint16_t mask, size_t x, size_t y, 
         uint64_t & sum, uint64_t & sumX, uint64_t & sumY) 
{ 
    int value = g_sum[mask]; 
    sum += value; 
    sumX += x * value + g_sumX[mask]; 
    sumY += y * value; 
} 

void GetCenter(const uint8_t * mask, size_t stride, size_t width, size_t height, 
    uint8_t index, double * centerX, double * centerY) 
{ 
    size_t alignedWidth = width & ~(sizeof(__m128i) - 1); 
    const __m128i _index = _mm_set1_epi8(index); 

    union PackedValue 
    { 
     uint8_t u8[4]; 
     uint16_t u16[2]; 
     uint32_t u32; 
    } _mask; 

    uint64_t sum = 0, sumX = 0, sumY = 0; 
    for(size_t y = 0; y < height; ++y) 
    { 
     size_t x = 0; 

     for(; x < alignedWidth; x += sizeof(__m128i)) 
     { 
      _mask.u32 = _mm_movemask_epi8(_mm_cmpeq_epi8(
         _mm_loadu_si128((__m128i*)(mask + x)), _index)); 
      AddMoments(_mask.u16[0], x, y, sum, sumX, sumY); 
     } 

     for(; x < width; ++x) 
     { 
      if(mask[x] == index) 
      { 
       sum++; 
       sumX += x; 
       sumY += y; 
      }    
     } 
     mask += stride; 
    } 
    *centerX = sum ? (double)sumX/sum : 0; 
    *centerY = sum ? (double)sumY/sum : 0; 
} 

Сравнение производительности для 1920х1080 изображения:

Как вы можете видеть выше, код с использованием 8-битных поисковых таблиц имеет более высокую производительность, чем код с использованием 16-разрядных таблиц поиска. Но в любом случае внешняя библиотека лучше, хотя и выполняет дополнительные вычисления моментов второго порядка.

Answer 2

Еще одна техника ускорения - кодирование по длине.

Вы можете разложить строки в горизонтальных прогонах, где активна маска. Вы можете обнаруживать пробеги «на лету» или предварительно компрометировать их и сохранять изображение в этой форме, если это имеет смысл.

Тогда прогон может быть накоплен в целом. Пусть начало бежать из (X, Y) и имеют длину L, а затем использовать

Sum+= L; 
SumX+= (2 * X + L + 1) * L; 
SumY+= Y * L; 

В конце концов, разделить SumX на 2.

Чем дольше прогоны, тем эффективнее трюк.

Используя SSE2 или более позднюю версию, попробуйте с инструкцией PMOVMSKB-Move Byte Mask.

Сначала сравните 16 пикселей маски с (реплицированным) значением индекса, чтобы получить 16 результатов сравнения. Затем упакуйте их в 16-битное число, используя магическую инструкцию.

Затем, используя две предварительно вычисленные таблицы поиска, выполните скопления в скалярном режиме.

Одна таблица поиска дает вам количество активных пикселей маски, а другая дает вам сумму активных пикселей маски, взвешенных по их абсциссе, то есть момент X.

Что-то вроде

int PackedValue= _mm_movemask_epi8(_mm_cmpeq_epi8(_mm_loadu_si128((__m128i*)&Mask[X]), ReplicatedIndex)); 
Sum+= Count[PackedValue]; 
SumX+= X * Count[PackedValue] + MomentX[PackedValue]; 
SumY+= Y * Count[PackedValue]; 

В зависимости от объема памяти, вы согласны потратить, таблицы подстановки могут иметь индексы байт (256 записей, используйте таблицу дважды) или индексы слово (65536 записей). В обоих случаях значения count и момента совпадают в одном байте (от 1 до 8/16 и от 0 до 28/120 соответственно).

AVX осуществление также возможен, упаковка 32 пикселей за раз. Однако таблицы поиска с индексами двойных слов кажутся необоснованными. :-)