在C语言中处理图像时算法的优化（图像卷积法

Question

我需要写一个程序，将图像读到一个数组中，进行卷积操作（锐化）并保存回文件（ppm）。我写了一个标准算法。

    unsigned char* imgBefore = malloc(height*(3*width)*sizeof(unsigned char));
    assert(fread(imgBefore, 3*width, height, inputFile) == height);
    unsigned char* imgAfter = malloc(height*(3*width)*sizeof(unsigned char));

    short ker[3][3] = {{0, -1, 0}, {-1, 5, -1}, {0, -1, 0}};

    unsigned char r, g, b;

    for(int y = 0; y < height; y++) {
        for(int x = 0; x < width; x++) {
            if(y == 0 || y == height-1 || x == 0 || x == width-1) {
                r = imgBefore[3*(width*y + x) + 0];
                g = imgBefore[3*(width*y + x) + 1];
                b = imgBefore[3*(width*y + x) + 2];
                imgAfter[3*(width*y + x) + 0]  = r;
                imgAfter[3*(width*y + x) + 1]  = g;
                imgAfter[3*(width*y + x) + 2]  = b;
                continue;
            }

            int rSum = 0, gSum = 0, bSum = 0, val;

            for(int dy = 0; dy < 3; dy++) { // kernels height
                int yy = 3*width*(y+dy-1);
                for(int dx = 0; dx < 3; dx++) { // kerenels width
                    int xx = 3*(x+dx-1);
                    val = ker[dy][dx];

                    rSum += val * imgBefore[yy + xx + 0];
                    gSum += val * imgBefore[yy + xx + 1];
                    bSum += val * imgBefore[yy + xx + 2];
                }
            }
            rSum = rSum < 0 ? 0 : (rSum > 255 ? 255 : rSum);
            gSum = gSum < 0 ? 0 : (gSum > 255 ? 255 : gSum);
            bSum = bSum < 0 ? 0 : (bSum > 255 ? 255 : bSum);

            imgAfter[3*(width*y + x) + 0] = rSum;
            imgAfter[3*(width*y + x) + 1] = gSum;
            imgAfter[3*(width*y + x) + 2] = bSum;
        }
    }

    fwrite(imgAfter, 3*width, height, outputFile);

接下来，我需要优化它与缓存的交互效果。在我看来，问题部分就在这段代码中。

for(int dy = 0; dy < 3; dy++) { // kernels height
    int yy = 3*width*(y+dy-1);
    for(int dx = 0; dx < 3; dx++) { // kerenels width
        int xx = 3*(x+dx-1);
        val = ker[dy][dx];

        rSum += val * imgBefore[yy + xx + 0];
        gSum += val * imgBefore[yy + xx + 1];
        bSum += val * imgBefore[yy + xx + 2];
    }
}

因为它首先加载了矩阵的一行，只用了3（9）个元素，然后转到下一行。这似乎与缓存的关系完全无效。

我该怎么做才能解决这个问题？

我也尝试过重用单个像素或整行。所有这些都只会使结果更加恶化（很有可能是我根本没有很好地实现像FIFO这样的结构，或者在错误的地方添加和读取了这些结构）。如果一个程序需要它，它应该是怎样的呢？评估时我使用：valgrind --tool=cachegrind --I1=32768,8,64 --D1=32768,8,64 --LL=1048576,16,64 .a.out。

如果有任何建议，我将不胜感激

Answer 1

对于大图像，你访问数据的方式对缓存不友好。这一点可以改进（例如使用分片）。

要知道，这个程序的计算时间不是受内存层次结构的限制，而是受数据布局效率低下导致的基于标量的慢序计算计算的限制。

下面是使用10次钢网函数调用的Valgrind分析结果。4608 x 3456 x 3 的形象。

--171871-- warning: L3 cache found, using its data for the LL simulation.
--171871-- warning: specified LL cache: line_size 64  assoc 12  total_size 9,437,184
--171871-- warning: simulated LL cache: line_size 64  assoc 18  total_size 9,437,184
==171871== 
==171871== I   refs:      30,437,416,378
==171871== I1  misses:               944
==171871== LLi misses:               938
==171871== I1  miss rate:           0.00%
==171871== LLi miss rate:           0.00%
==171871== 
==171871== D   refs:       7,526,412,742  (7,000,797,573 rd   + 525,615,169 wr)
==171871== D1  misses:        30,599,898  (   22,387,768 rd   +   8,212,130 wr)
==171871== LLd misses:        15,679,220  (    7,467,138 rd   +   8,212,082 wr)
==171871== D1  miss rate:            0.4% (          0.3%     +         1.6%  )
==171871== LLd miss rate:            0.2% (          0.1%     +         1.6%  )
==171871== 
==171871== LL refs:           30,600,842  (   22,388,712 rd   +   8,212,130 wr)
==171871== LL misses:         15,680,158  (    7,468,076 rd   +   8,212,082 wr)
==171871== LL miss rate:             0.0% (          0.0%     +         1.6%  )

首先，我们可以看到，有很多的。D1 cache reference 每写一个像素47个，由于标量访问）。

的数量。D1 misses 和 D1 miss rate 似乎非常小，所以我们可以认为缓存访问是好的。但事实并非如此。Valgrind的结果非常具有误导性，因为 D1 cache reference 工作在1字节的标量粒度上，而 D1 misses 工作在64字节的缓存行粒度。

如果我们更仔细地看分析，可以看到几乎所有的写入数据都会造成D1缓存的错过。而且，从L1读取的数据中，有20%需要从LLC中（重新）加载（对于从D1缓存读取的7GB，有1.43GB是从LLC中加载的）。

问题来自于图像的大行。事实上，1行图像需要 4608 * 3 = 13824 字节，写一行需要3行。因此，从D1缓存读取41472个字节，要写13824个字节。理论上，如果D1缓存足够大，应该可以从D1缓存中重用27648个。但是，D1缓存太小，无法容纳所有的读写数据（大致是因为 41472+13824 > 32768).

我们可以用两种方法来改善这种情况。

• 平铺直叙：其思路是将 x 循环绑定在一个小的固定大小（例如1024），这样它将计算图像的垂直块，然后添加一个新的包围的图像。xBlock 循环（其中包括 y 和 x 循环），在垂直区块上迭代。这种方法的好处是提高了D1缓存的重用性，但在现代处理器上仍可能造成一些不必要的缓存遗漏（因为加载的行已经不连续了）。您可以通过在下面的数据块上执行平铺来缓解这种情况。y 维度或通过调整数据布局。下面是一个例子。

const int blockSize = 1024;
for(int xBlock = 0; xBlock < width; xBlock+=blockSize) {
    for(int y = 0; y < height; y++) {
        int xMax = xBlock + blockSize;
        if(xMax > width)
            xMax = width;
        for(int x = xBlock; x < xMax; x++) {
            // Unchanged part of the code
        }
    }
}

• 非时态存储：向缓存层次结构写入数据会造成缓存污染，导致额外的缓存遗漏（因为写入数据需要驱逐之前读取的缓存行）。它还可能会在某些架构上造成更多的负载（比如像Intel处理器上）。我们可以告诉处理器直接将数据写入内存。一些编译器提供了 #pragma 指令来实现这个功能（在你的代码中手动操作有点麻烦）。只有当图像太大，无法放入LLC缓存中或以后无法重用时，才可以这样做。请看 此处 以获取更多信息。

下面是应用平铺法后的Valgrind结果。

==183432== D   refs:       7,527,450,132  (7,001,731,093 rd   + 525,719,039 wr)
==183432== D1  misses:        16,025,288  (    7,640,368 rd   +   8,384,920 wr)
==183432== LLd misses:        16,024,790  (    7,639,918 rd   +   8,384,872 wr)

我们可以看到，现在D1缓存的失误减少了2倍（读取次数减少了3倍）。