利用 CPU 流水线加快数据处理

最近看 OpenCV 源码时注意到一个有意思的地方：

template<typename T, size_t BinsOnStack = 0u>
static double getThreshVal_Otsu( const Mat& _src, const Size& size)
{
    const int N = std::numeric_limits<T>::max() + 1;
    int i, j;
    #if CV_ENABLE_UNROLLED
    AutoBuffer<int, 4 * BinsOnStack> hBuf(4 * N); // 分配了4倍 buffer
    #else
    AutoBuffer<int, BinsOnStack> hBuf(N);
    #endif
    memset(hBuf.data(), 0, hBuf.size() * sizeof(int));
    int* h = hBuf.data();
    #if CV_ENABLE_UNROLLED
    int* h_unrolled[3] = {h + N, h + 2 * N, h + 3 * N };
    #endif
    for( i = 0; i < size.height; i++ )
    { // 统计灰度直方图
        const T* src = _src.ptr<T>(i, 0);
        j = 0;
        #if CV_ENABLE_UNROLLED
        for( ; j <= size.width - 4; j += 4 )
        { // 行方向上的迭代次数减少为 1/4
            int v0 = src[j], v1 = src[j+1];
            h[v0]++; h_unrolled[0][v1]++;
            v0 = src[j+2]; v1 = src[j+3];
            h_unrolled[1][v0]++; h_unrolled[2][v1]++;
        }
        #endif
        for( ; j < size.width; j++ )
            h[src[j]]++;
    }

    double mu = 0, scale = 1./(size.width*size.height);
    for( i = 0; i < N; i++ )
    {
        #if CV_ENABLE_UNROLLED
        h[i] += h_unrolled[0][i] + h_unrolled[1][i] + h_unrolled[2][i];
        #endif
        mu += i*(double)h[i];
    }

    mu *= scale;
    double mu1 = 0, q1 = 0;
    double max_sigma = 0, max_val = 0;

    for(i = 0; i < N; i++ )
    {
        double p_i, q2, mu2, sigma;

        p_i = h[i]*scale;
        mu1 *= q1;
        q1 += p_i;
        q2 = 1. - q1;

        if( std::min(q1,q2) < FLT_EPSILON || std::max(q1,q2) > 1. - FLT_EPSILON )
            continue;

        mu1 = (mu1 + i*p_i)/q1;
        mu2 = (mu - q1*mu1)/q2;
        sigma = q1*q2*(mu1 - mu2)*(mu1 - mu2);
        if( sigma > max_sigma )
        {
            max_sigma = sigma;
            max_val = i;
        }
    }
    return max_val;
}

这是大津法求二值化阈值的过程，里面有个CV_ENABLE_UNROLLED宏（默认为 1，即开启）。经过测试，我发现开启该宏有比较明显的加速效果。当这个宏开启时，行方向上的迭代次数减少为 1/4，加上宏名中带有“UNROLLED”，起初我以为这里应该是利用 loop unrolling（循环展开）或者 loop tiling（循环分块）达到加速效果的吧，但分配 4 倍的 buffer 有必要吗？

于是我去掉了 4 倍 buffer，将循环部分改成：

for (; j <= image.cols - 4; j += 4) {
    h[src[j]]++;
    h[src[j + 1]]++;
    h[src[j + 2]]++;
    h[src[j + 3]]++;
}

神奇的是，加速效果完全消失。

百思不得其解，网上也很难找到有用的资料。后来跟同事探讨的过程中，经提醒还有个 CPU 流水线的概念 (图片来自 Eigen CGLibs 2013. Giugno Pisa. P37)：

CPU 流水线（Pipelining）大概的意思就是：CPU 使用不同的处理单元进行不同的数据操作（比如从内存取数据、数据相加、存数据到内存），对同一内存单元的不同操作只能串行执行，但对不同内存单元的不同操作可以自动并行！

回过头看 otsu 算法，由于图像中相邻像素很可能有相同的像素值，即 h[src[j + ?]] 很可能指向同一个内存单元，这种情况下循环中前后两条指令就无法利用流水线，只能串行执行。但如果使用 4 个 buffer，则可以确保循环中 4 条语句使用不同的内存单元，前一条语句进行加法操作时，后一条语句可以同时进行数据读取操作。

速度测试程序如下：

#include <chrono>
#include <opencv2/opencv.hpp>

using namespace cv;
using namespace std;
using namespace chrono;

template <bool ENABLE_UNROLL>
static double pixCount(const Mat& image)
{
    constexpr int N = 256;
    constexpr int buff_size = ENABLE_UNROLL ? N * 4 : N;
    int i, j;
    std::array<int, buff_size> hBuf;
    memset(hBuf.data(), 0, hBuf.size() * sizeof(int));
    int* h = hBuf.data();
    int* h_unrolled[3] = { h + N, h + 2 * N, h + 3 * N };
    for (i = 0; i < image.rows; i++) {
        const uchar* src = image.ptr<uchar>(i, 0);
        j = 0;
        for (; j <= image.cols - 4; j += 4) {
            h[src[j]]++;
            if constexpr (ENABLE_UNROLL) {
                h_unrolled[0][src[j + 1]]++;
                h_unrolled[1][src[j + 2]]++;
                h_unrolled[2][src[j + 3]]++;
            } else {
                h[src[j + 1]]++;
                h[src[j + 2]]++;
                h[src[j + 3]]++;
            }
        }
        for (; j < image.cols; j++)
            h[src[j]]++;
    }

    double mu = 0;
    for (i = 0; i < N; i++) {
        if constexpr (ENABLE_UNROLL) {
            h[i] += h_unrolled[0][i] + h_unrolled[1][i] + h_unrolled[2][i];
        }
        mu += i * (double)h[i];
    }
    return mu;
}

int main(int argc, char** argv)
{
    Mat image = imread(argv[1], IMREAD_GRAYSCALE);
    cout << "image size: " << image.cols << "x" << image.rows << endl;
    auto time_start = system_clock::now();
    auto res = pixCount<false>(image);
    auto time_cost = duration_cast<nanoseconds>(system_clock::now() - time_start).count();
    cout << "no unrolling: " << time_cost * 1e-6 << " ms \tresult: " << res << endl;
    time_start = system_clock::now();
    res = pixCount<true>(image);
    time_cost = duration_cast<nanoseconds>(system_clock::now() - time_start).count();
    cout << "enable unrolling: " << time_cost * 1e-6 << " ms \tresult: " << res << endl;
}

执行结果：

image size: 3200x2400
no unrolling: 9.12246 ms 	    result: 1.16491e+08
enable unrolling: 3.35144 ms 	result: 1.16491e+08

当图像中相邻像素值相同的概率很大时（比如存在纯色背景或者图像内容简单），利用 CPU 流水线的处理方式加速效果非常显著。