图像矩的快速计算

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

有一种方法可以大大（超过十倍）提高该算法的性能。 为此，您需要使用 CPU 的 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;
}

附：使用外部库 (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;
}

1920x1080 图像的性能比较：

Base version: 8.261 ms; 
1-st optimization:0.363 ms (in 22 times faster);
2-nd optimization:0.280 ms (in 29 times faster);
3-rd optimization:0.299 ms (in 27 times faster);
4-th optimization:0.325 ms (in 25 times faster);

正如您在上面看到的那样，使用 8 位查找表的代码比使用 16 位查找表的代码具有更好的性能。但无论如何，外部库更好，尽管它执行二阶矩的额外计算。

Answer 2

另一种加速技术是游程编码。

您可以在掩码处于活动状态的水平运行中分解行。如果有意义的话，您可以动态检测运行，或预先计算它们并以该形式存储图像。

然后一个运行可以作为一个整体进行累积。让运行从(X, Y) 开始，长度为L，然后使用

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

最后，将SumX 除以2。

运行时间越长，技巧越有效。

---

使用 SSE2 或更高版本，您可以尝试使用指令 PMOVMSKB—移动字节掩码。

首先将 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 个条目）。在这两种情况下，计数和矩值都适合单个字节（分别为 1 到 8/16 和 0 到 28/120）。

AVX 实现也是可能的，一次打包 32 个像素。但是，带有双字索引的查找表似乎不合理。 :-)