在 C 中处理图像时优化算法(图像卷积)

Optimization of the algorithm when working with images in C (image convolution)

我需要编写一个程序,将图像读入数组,执行卷积运算(锐化)并将其保存回文件 (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

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

您访问数据的方式对于大图像不是缓存友好的。这可以改进(例如使用平铺)。

请注意,此程序的计算时间不受内存层次结构的限制,而是由低效的数据布局导致的基于顺序标量的缓慢计算。

这是使用 4608 x 3456 x 3 图片对 10 次模板函数调用的 Valgrind 分析结果:

--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 missesD1 miss rate的数量似乎很少,所以我们可以认为缓存访问很好。但事实并非如此。 Valgrind 结果非常具有误导性,因为 D1 cache reference 以 1 字节标量粒度工作,而 D1 misses 以 64 字节缓存行粒度工作。

如果我们仔细分析一下,我们可以看到几乎所有写入的数据都会导致D1缓存未命中。此外,从 L1 读取的数据中有 20% 需要从 LLC(重新)加载(对于从 D1 缓存读取的 7 GB,从 LLC 加载 1.43 GB)。

问题来自图像的大线条。事实上,1 条图像行需要 4608 * 3 = 13824 字节,写一行需要 3 行。因此,从 D1 缓存中读取 41472 字节写入 13824 字节。理论上,如果 D1 缓存足够大,应该从 D1 缓存中重用 27648。但是D1缓存太小,放不下所有read/written个数据(大概是因为41472+13824 > 32768)。

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

  • 平铺:想法是将 x 循环边界设置为较小的固定大小(例如 1024),以便计算图像的垂直块,然后添加新的包含 xBlock 循环(包括 yx 循环)迭代垂直块。这种方法具有提高 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 缓存或以后不再使用时,这才有用。有关详细信息,请参阅 here

这是应用平铺方法后的 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 倍)。