在 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 misses
和D1 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
循环(包括 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 缓存或以后不再使用时,这才有用。有关详细信息,请参阅 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 倍)。
我需要编写一个程序,将图像读入数组,执行卷积运算(锐化)并将其保存回文件 (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 misses
和D1 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
循环(包括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 缓存或以后不再使用时,这才有用。有关详细信息,请参阅 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 倍)。