CSAPP LAB 实验

实验介绍

图像处理中存在很多函数，可以对这些函数进行优化。本实验主要关注两种图像处理操作：

  • 旋转：对图像逆时针旋转 90 度
  • 平滑：对图像进行模糊操作

图像用二维矩阵 M 表示，MijM_{ij} 表示图像 M 的第 (i,j)(i,j) 像素的值，像素值用红，绿，蓝表示。

旋转操作用下面 2 个操作表示：

  • Transpose：对第（i,j）个像素对，MijM_{ij} 和 MjiM_{ji} 交换
  • Exchange rows：行 i 和行 N-1-i 交换

平滑操作：每个像素用周围像素值的平均值表示（相当于进行均值滤波）。比如：

M2[1][1]=∑i=02∑j=02M1[i][j]9M2[N−1][N−1]=∑i=N−2N−1∑j=N−2N−1M1[i][j]4\begin{aligned} M_2[1][1] = & \frac{\sum_{i=0}^2 \sum_{j=0}^2 M_1[i][j]}{9} \\ M_2[N-1][N-1] = & \frac{\sum_{i=N-2}^{N-1} \sum_{j=N-2}^{N-1} M_1[i][j]}{4} \end{aligned}

image.png

编码规则：

  • 只能用 ANSI C，不能用嵌入式汇编
  • 不能修改测量时间的机制（CPE）
  • 只能修改 kernels.c，可以定义宏，全局变量，函数

实验开始

知识回顾

回顾一下 CSAPP 第 5 第 6 章学到了哪些优化代码的小技巧：

  • 消除冗余的函数调用。比如避免在 for 循环里用 strlen。
  • 消除不必要的内存引用。比如引入临时变量来把中间结果保存到寄存器里，在全部计算完成后把最终结果存到数组或全局变量里。
  • 循环展开，降低判断语句和增减循环变量的开销。
  • 累积变量和重新组合，提高指令并行性。
  • 功能性风格重写条件操作，即用三元运算符。
  • 提高空间局部性，尽量按照数组在内存里存储的顺序，以 1 为步长进行读取。
  • 提高时间局部性，一旦在内存里读出一个变量，就尽可能频繁且集中的使用它。

Cache 的理论依据：

  • 时间局部性：如果程序中的某条指令一旦执行，不久以后该指令可能再次执行；如果某数据被访问过，不久以后该数据可能再次被访问。也就是说：被引用过一次的存储器位置在未来会被多次引用，比如循环操作。
  • 空间局部性：一旦程序访问了某个存储单元，在不久之后，其附近的存储单元也将被访问。也就是说：如果一个存储器的位置被引用，那么将来他附近的位置也会被引用，比如 cache。
  • Cache 有读命中和写命中，写不命中开销比读不命中大。因为在计算机缓存机制中，写操作通常比读操作更为复杂和成本较高。

准备工作

解压作业包：tar xvf perflab-handout.tar。

只能允许修改 kernels.c 文件。我们首先在 kernels.c 填上自己的信息，程序才能成功启动：

1
2
3
4
5
6
7
8
9
10
11
12
/*
 * Please fill in the following team struct
 */
team_t team = {
    "casa", /* Team name */

    "",                /* First member full name */
    "", /* First member email address */

    "", /* Second member full name (leave blank if none) */
    ""  /* Second member email addr (leave blank if none) */
};

测试和运行：

1
2
make driver
./driver

得到的其中一条示例结果为：

1
2
3
4
5
Smooth: Version = naive_smooth: Naive baseline implementation:
Dim             32      64      128     256     512     Mean
Your CPEs       27.0    28.0    26.3    29.2    24.1
Baseline CPEs   695.0   698.0   702.0   717.0   722.0
Speedup         25.7    25.0    26.6    24.6    29.9    26.3

有时候运行会抽风，计算的 CPE 出现负数。并且提示：Fatal Error: Non-positive CPE value… 这时候重新运行一遍就行。另外每次运行 driver 得到的分数都会不一样。

第一行表示执行的方法与函数信息，我们的目标让自己的函数的 Your CPEs 尽可能低，让 Speedup 的 Mean（平均值）比使用默认方法要高。

CPE or Cycles per Element。如果一个函数使用 C 个周期去执行一个大小为 N*N 的图像，那么 CPE=C/(N*N),因此 CPE 越小越好。

在 kernel.c 中有 rotate 的实现方法，我们需要添加自己的 rotate 实现，函数名自己命名，然后通过以下方式注册函数：

1
2
3
4
void register_rotate_functions() { 
	add_rotate_function(&rotate, rotate_descr); 
	add_rotate_function(&my_rotate，my_rotate_descr);
}

每个像素的结构：

1
2
3
4
5
typedef struct { 
	unsigned short red; /* R value */ 
	unsigned short green; /* G value */ 
	unsigned short blue; /* B value */ 
} pixel

图像用一维数组表示，第（i,j）个像素表示为 I[RIDX(i,j,n)]，n 为图像的维数。

1
#define RIDX(i,j,n) ((i)*(n)+(j))

图像旋转

原函数的实现为：

1
2
3
4
5
6
7
8
void naive_rotate(int dim, pixel *src, pixel *dst)
{
    int i, j;

    for (i = 0; i < dim; i++)
        for (j = 0; j < dim; j++)
            dst[RIDX(dim - 1 - j, i, dim)] = src[RIDX(i, j, dim)];
}

减少写数据的迭代步长

利用知识点：写不命中开销比读不命中大。我们应该优先对写入像素点的索引进行优化。

观察原函数得知，dst 在进行写时，单次迭代的步长很大。以下是针对这点进行改进的三个版本，它们都比原始方法的性能更好。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
/*
	单纯交换两个for循环的i和j
*/
char rotateij_rev_descr[] = "rotate:  单纯 i,j 调换(遍历顺序调换)";
void rotateij_rev(int dim, pixel *src, pixel *dst)
{
    int i, j;

    for (j = 0; j < dim; j++)
        for (i = 0; i < dim; i++)
            dst[RIDX(dim - 1 - j, i, dim)] = src[RIDX(i, j, dim)];
}


/*
	改写取dst索引的方式，按顺序写入dst。
	索引的推导方法可以通过自己手动模拟，写出src对应的索引式子：jn+n-1-i = (j+1)n-1-i
*/
char rotateij_descr[] = "rotate:  单纯 i,j 调换";
void rotateij(int dim, pixel *src, pixel *dst)
{
    int i, j;

    for (i = 0; i < dim; i++)
        for (j = 0; j < dim; j++)
            dst[RIDX(i, j, dim)] = src[RIDX(j+1, -1-i, dim)];
}

/*
	在rotateij的基础上把索引变量提取出来，尽量不使用乘法。使用指针代替RIDX进行数组访问。
*/
char rotateij_local_descr[] = "rotate: i,j 调换 + 使用局部变量 ";
void rotateij_local(int dim, pixel *src, pixel *dst)
{
    int i, j;
    int didx = 0;
    int osidx = RIDX(1, -1, dim);
    int sidx = osidx--;
    for (i = 0; i < dim; i++)
    {
        for (j = 0; j < dim; j++)
        {
            dst[didx++] = src[sidx];
            sidx += dim;
        }
        sidx = osidx--;
    }
}

结果如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
Rotate: Version = rotate:  单纯 i,j 调换(遍历顺序调换):
Dim             64      128     256     512     1024    Mean
Your CPEs       1.0     1.1     2.0     2.3     5.8
Baseline CPEs   14.7    40.1    46.4    65.9    94.5
Speedup         14.3    35.2    23.3    28.6    16.4    22.3

Rotate: Version = rotate:  单纯 i,j 调换:
Dim             64      128     256     512     1024    Mean
Your CPEs       1.3     1.2     1.9     2.6     5.8
Baseline CPEs   14.7    40.1    46.4    65.9    94.5
Speedup         11.2    32.1    24.3    25.1    16.2    20.4

Rotate: Version = rotate: i,j 调换 + 使用局部变量 :
Dim             64      128     256     512     1024    Mean
Your CPEs       1.3     1.2     1.8     2.5     5.6
Baseline CPEs   14.7    40.1    46.4    65.9    94.5
Speedup         11.3    32.9    25.6    26.4    16.9    21.2

观察可知 rotateij_rev 效果更好。由于历史原因（由于我懒得改代码了），我这里使用了 rotateij_local 方法。读者可随意选取方法，不影响后续的改进。

矩阵分块计算

观察到测试的图像长宽都是 16 或 32 的倍数。我们可以考虑矩阵分块。

矩阵分块可以减少缓存的不命中率。空间局部性更好。这里测试了 32*32 分块和 16*16 分块代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
char rotate_blc32_descr[] = "rotate: 32 大小的块";
void rotate_blc32(int dim, pixel *src, pixel *dst)
{
    // naive_rotate(dim, src, dst);
    int i, j, I, J;
    int block_size = 32;
    for (I = 0; I < dim; I += block_size)
    {
        for (J = 0; J < dim; J += block_size)
        {
            for (i = 0; i < block_size; i++)
            {
                for (j = 0; j < block_size; j++)
                {
                    dst[RIDX(i + I, j + J, dim)] = src[RIDX(J + j + 1, -1 - i - I, dim)];
                }
            }
        }
    }
}

char rotate_blc16_descr[] = "rotate: 16 大小的块";
void rotate_blc16(int dim, pixel *src, pixel *dst)
{
    // naive_rotate(dim, src, dst);
    int i, j, I, J;
    int block_size = 16;
    for (I = 0; I < dim; I += block_size)
    {
        for (J = 0; J < dim; J += block_size)
        {
            for (i = 0; i < block_size; i++)
            {
                for (j = 0; j < block_size; j++)
                {
                    dst[RIDX(i + I, j + J, dim)] = src[RIDX(J + j + 1, -1 - i - I, dim)];
                }
            }
        }
    }
}

多次测试表明按 16*16 分块效果会好一点。

1
2
3
4
5
6
7
8
9
10
11
Rotate: Version = rotate: 32 大小的块:
Dim             64      128     256     512     1024    Mean
Your CPEs       1.6     1.4     1.5     1.7     2.5
Baseline CPEs   14.7    40.1    46.4    65.9    94.5
Speedup         9.0     28.1    31.4    38.6    37.5    25.8

Rotate: Version = rotate: 16 大小的块:
Dim             64      128     256     512     1024    Mean
Your CPEs       1.4     1.5     1.4     1.8     2.5
Baseline CPEs   14.7    40.1    46.4    65.9    94.5
Speedup         10.2    27.5    32.7    36.0    38.0    26.3

关于矩阵分块的具体原理推荐阅读：

循环展开

循环展开是一种牺牲程序的尺寸来加快程序的执行速度的优化方法。可以由程序员完成，也可由编译器自动优化完成。循环展开的本质是，利用 CPU 指令级并行，来降低循环的开销，当然，同时也有利于指令流水线的高效调度。

优点：

  • 提高缓存命中（cache hit）率，增加循环体内语句并发执行的可能性（需要循环体内语句不相关）；
  • 减少分支预测失败的可能性，提高性能

缺点：

  • 程序代码膨胀、代码可读性降低
  • 消耗较多寄存器缓存

在分块的基础上进行循环展开：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
char rotate_blc16_16_descr[] = "rotate: 16 块 16*1 循环展开";
void rotate_blc16_16(int dim, pixel *src, pixel *dst)
{
    int i, I, J;
    int block_size = 16;
    int didx, sidx;
    for (I = 0; I < dim; I += block_size)
    {
        for (J = 0; J < dim; J += block_size)
        {
            for (i = 0; i < block_size; i++)
            {
                didx = RIDX(i + I, J, dim);
                sidx = RIDX(J + 1, -1 - i - I, dim);
                dst[didx++] = src[sidx];
                sidx += dim;
                dst[didx++] = src[sidx];
                sidx += dim;
                dst[didx++] = src[sidx];
                sidx += dim;
                dst[didx++] = src[sidx];
                sidx += dim;
                dst[didx++] = src[sidx];
                sidx += dim;
                dst[didx++] = src[sidx];
                sidx += dim;
                dst[didx++] = src[sidx];
                sidx += dim;
                dst[didx++] = src[sidx];
                sidx += dim;
                dst[didx++] = src[sidx];
                sidx += dim;
                dst[didx++] = src[sidx];
                sidx += dim;
                dst[didx++] = src[sidx];
                sidx += dim;
                dst[didx++] = src[sidx];
                sidx += dim;
                dst[didx++] = src[sidx];
                sidx += dim;
                dst[didx++] = src[sidx];
                sidx += dim;
                dst[didx++] = src[sidx];
                sidx += dim;
                dst[didx++] = src[sidx];
                sidx += dim;
            }
        }
    }
}
1
2
3
4
5
Rotate: Version = rotate: 16 块 16*1 循环展开:
Dim             64      128     256     512     1024    Mean
Your CPEs       1.2     1.2     1.3     1.5     2.3
Baseline CPEs   14.7    40.1    46.4    65.9    94.5
Speedup         12.4    33.8    34.6    45.3    40.4    30.5

平滑操作

原函数存在的问题是频繁的调用 avg 函数、min、max 等各种函数开销较大。

我们可以考虑循环展开或者消除函数调用，进一步，可以将重复利用的数据存储在了数组之中。

思路：

  • 不使用其他任何函数，自己实现平滑计算
  • 对四个点、四条边和中间内容分别处理
  • 充分利用统计先前的统计结果统计像素值的总和

temp_pixel_group 用于统计列的像素之和。为了便于迭代操作和空间利用，temp_pixel_group 安排方式为：

image.png

循环展开后的代码看着眼花，调试 Bug 更是折磨人。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
// 定义自己的数据结构。注意成员的类型是int类型。
typedef struct
{
    int red;
    int green;
    int blue;
} my_pixel_sum;

char smooth_my_descr[] = "smooth: 我的实现版本";
void smooth_my(int dim, pixel *src, pixel *dst){
    int i, j;

    // 四点处理
    int q1, q2, q3, q4, idx;
    q1 = 0;
    q2 = 1;
    q3 = dim;
    q4 = dim + 1;
    idx = 0;
    dst[idx].blue = (src[q1].blue + src[q2].blue + src[q3].blue + src[q4].blue) / 4;
    dst[idx].green = (src[q1].green + src[q2].green + src[q3].green + src[q4].green) / 4;
    dst[idx].red = (src[q1].red + src[q2].red + src[q3].red + src[q4].red) / 4;

    q1 = dim - 2;
    q2 = dim - 1;
    q3 = dim + dim - 2;
    q4 = dim + dim - 1;
    idx = dim - 1;
    dst[idx].blue = (src[q1].blue + src[q2].blue + src[q3].blue + src[q4].blue) / 4;
    dst[idx].green = (src[q1].green + src[q2].green + src[q3].green + src[q4].green) / 4;
    dst[idx].red = (src[q1].red + src[q2].red + src[q3].red + src[q4].red) / 4;

    q1 = (dim - 2) * dim;
    q2 = (dim - 2) * dim + 1;
    q3 = (dim - 1) * dim;
    q4 = (dim - 1) * dim + 1;
    idx = (dim - 1) * dim;
    dst[idx].blue = (src[q1].blue + src[q2].blue + src[q3].blue + src[q4].blue) / 4;
    dst[idx].green = (src[q1].green + src[q2].green + src[q3].green + src[q4].green) / 4;
    dst[idx].red = (src[q1].red + src[q2].red + src[q3].red + src[q4].red) / 4;

    q1 = dim * dim - 2 - dim;
    q2 = dim * dim - 1 - dim;
    q3 = dim * dim - 2;
    q4 = dim * dim - 1;
    idx = dim * dim - 1;
    dst[idx].blue = (src[q1].blue + src[q2].blue + src[q3].blue + src[q4].blue) / 4;
    dst[idx].green = (src[q1].green + src[q2].green + src[q3].green + src[q4].green) / 4;
    dst[idx].red = (src[q1].red + src[q2].red + src[q3].red + src[q4].red) / 4;

    my_pixel_sum *temp_pixel_group = (my_pixel_sum *)malloc(sizeof(my_pixel_sum) * 3);

    // 上边处理
    // 初始化temp_pixel_group
    for (i = 0; i <= 1; i++)
    {
        temp_pixel_group[i].blue = src[RIDX(0, i, dim)].blue + src[RIDX(1, i, dim)].blue;
        temp_pixel_group[i].green = src[RIDX(0, i, dim)].green + src[RIDX(1, i, dim)].green;
        temp_pixel_group[i].red = src[RIDX(0, i, dim)].red + src[RIDX(1, i, dim)].red;
    }
    for (j = 1; j <= dim - 2; j++)
    {
        q1 = j + 1;
        temp_pixel_group[q1 % 3].blue = src[RIDX(0, q1, dim)].blue + src[RIDX(1, q1, dim)].blue;
        temp_pixel_group[q1 % 3].green = src[RIDX(0, q1, dim)].green + src[RIDX(1, q1, dim)].green;
        temp_pixel_group[q1 % 3].red = src[RIDX(0, q1, dim)].red + src[RIDX(1, q1, dim)].red;
        dst[RIDX(0, j, dim)].blue = (temp_pixel_group[0].blue + temp_pixel_group[1].blue + temp_pixel_group[2].blue) / 6;
        dst[RIDX(0, j, dim)].red = (temp_pixel_group[0].red + temp_pixel_group[1].red + temp_pixel_group[2].red) / 6;
        dst[RIDX(0, j, dim)].green = (temp_pixel_group[0].green + temp_pixel_group[1].green + temp_pixel_group[2].green) / 6;
    }

    // 下边处理
    // 初始化temp_pixel_group
    for (i = 0; i <= 1; i++)
    {
        temp_pixel_group[i].blue = src[RIDX(dim - 1, i, dim)].blue + src[RIDX(dim - 2, i, dim)].blue;
        temp_pixel_group[i].green = src[RIDX(dim - 1, i, dim)].green + src[RIDX(dim - 2, i, dim)].green;
        temp_pixel_group[i].red = src[RIDX(dim - 1, i, dim)].red + src[RIDX(dim - 2, i, dim)].red;
    }
    for (j = 1; j <= dim - 2; j++)
    {
        q1 = j + 1;
        // q2 = q1 + dim;
        temp_pixel_group[q1 % 3].blue = src[RIDX(dim - 1, q1, dim)].blue + src[RIDX(dim - 2, q1, dim)].blue;
        temp_pixel_group[q1 % 3].green = src[RIDX(dim - 1, q1, dim)].green + src[RIDX(dim - 2, q1, dim)].green;
        temp_pixel_group[q1 % 3].red = src[RIDX(dim - 1, q1, dim)].red + src[RIDX(dim - 2, q1, dim)].red;
        dst[RIDX(dim - 1, j, dim)].blue = (temp_pixel_group[0].blue + temp_pixel_group[1].blue + temp_pixel_group[2].blue) / 6;
        dst[RIDX(dim - 1, j, dim)].red = (temp_pixel_group[0].red + temp_pixel_group[1].red + temp_pixel_group[2].red) / 6;
        dst[RIDX(dim - 1, j, dim)].green = (temp_pixel_group[0].green + temp_pixel_group[1].green + temp_pixel_group[2].green) / 6;
    }

    // 左边处理
    // 初始化temp_pixel_group
    for (i = 0; i <= 1; i++)
    {
	    
        temp_pixel_group[i].blue = src[RIDX(i, 0, dim)].blue + src[RIDX(i, 1, dim)].blue;
        temp_pixel_group[i].green = src[RIDX(i, 0, dim)].green + src[RIDX(i, 1, dim)].green;
        temp_pixel_group[i].red = src[RIDX(i, 0, dim)].red + src[RIDX(i, 1, dim)].red;
    }
    for (i = 1; i <= dim - 2; i++)
    {
        q1 = i + 1;
        // q2 = q1 + dim;
        temp_pixel_group[q1 % 3].blue = src[RIDX(q1, 0, dim)].blue + src[RIDX(q1, 1, dim)].blue;
        temp_pixel_group[q1 % 3].green = src[RIDX(q1, 0, dim)].green + src[RIDX(q1, 1, dim)].green;
        temp_pixel_group[q1 % 3].red = src[RIDX(q1, 0, dim)].red + src[RIDX(q1, 1, dim)].red;
        dst[RIDX(i, 0, dim)].blue = (temp_pixel_group[0].blue + temp_pixel_group[1].blue + temp_pixel_group[2].blue) / 6;
        dst[RIDX(i, 0, dim)].red = (temp_pixel_group[0].red + temp_pixel_group[1].red + temp_pixel_group[2].red) / 6;
        dst[RIDX(i, 0, dim)].green = (temp_pixel_group[0].green + temp_pixel_group[1].green + temp_pixel_group[2].green) / 6;
    }

    // 右边处理
    // 初始化temp_pixel_group
    for (i = 0; i <= 1; i++)
    {
        temp_pixel_group[i].blue = src[RIDX(i, dim - 1, dim)].blue + src[RIDX(i, dim - 2, dim)].blue;
        temp_pixel_group[i].green = src[RIDX(i, dim - 1, dim)].green + src[RIDX(i, dim - 2, dim)].green;
        temp_pixel_group[i].red = src[RIDX(i, dim - 1, dim)].red + src[RIDX(i, dim - 2, dim)].red;
    }
    for (i = 1; i <= dim - 2; i++)
    {
        q1 = i + 1;
        temp_pixel_group[q1 % 3].blue = src[RIDX(q1, dim - 1, dim)].blue + src[RIDX(q1, dim - 2, dim)].blue;
        temp_pixel_group[q1 % 3].green = src[RIDX(q1, dim - 1, dim)].green + src[RIDX(q1, dim - 2, dim)].green;
        temp_pixel_group[q1 % 3].red = src[RIDX(q1, dim - 1, dim)].red + src[RIDX(q1, dim - 2, dim)].red;
        dst[RIDX(i, dim - 1, dim)].blue = (temp_pixel_group[0].blue + temp_pixel_group[1].blue + temp_pixel_group[2].blue) / 6;
        dst[RIDX(i, dim - 1, dim)].red = (temp_pixel_group[0].red + temp_pixel_group[1].red + temp_pixel_group[2].red) / 6;
        dst[RIDX(i, dim - 1, dim)].green = (temp_pixel_group[0].green + temp_pixel_group[1].green + temp_pixel_group[2].green) / 6;
    }

    // 中间通用
    int k;
    for (i = 1; i <= dim - 2; i++)
    {
        // 换行初始化
        for (k = 0; k <= 2; k++)
        {
            temp_pixel_group[k].blue = src[RIDX(i-1, k, dim)].blue + src[RIDX(i, k, dim)].blue + src[RIDX(i+1, k, dim)].blue;
            temp_pixel_group[k].green = src[RIDX(i-1, k, dim)].green + src[RIDX(i, k, dim)].green + src[RIDX(i+1, k, dim)].green;
            temp_pixel_group[k].red = src[RIDX(i-1, k, dim)].red + src[RIDX(i, k, dim)].red + src[RIDX(i+1, k, dim)].red;
        }
        for (j = 1; j <= dim - 2; j++)
        {
            q1 = j+1;
            temp_pixel_group[q1 % 3].blue = src[RIDX(i-1, j+1, dim)].blue + src[RIDX(i, j+1, dim)].blue+ src[RIDX(i+1, j+1, dim)].blue;
            temp_pixel_group[q1 % 3].green = src[RIDX(i-1, j+1, dim)].green + src[RIDX(i, j+1, dim)].green+ src[RIDX(i+1, j+1, dim)].green;
            temp_pixel_group[q1 % 3].red = src[RIDX(i-1, j+1, dim)].red + src[RIDX(i, j+1, dim)].red+ src[RIDX(i+1, j+1, dim)].red;
            dst[RIDX(i, j, dim)].blue = (temp_pixel_group[0].blue + temp_pixel_group[1].blue + temp_pixel_group[2].blue) / 9;
            dst[RIDX(i, j, dim)].red = (temp_pixel_group[0].red + temp_pixel_group[1].red + temp_pixel_group[2].red) / 9;
            dst[RIDX(i, j, dim)].green = (temp_pixel_group[0].green + temp_pixel_group[1].green + temp_pixel_group[2].green) / 9;
        }
    }
}

后续可进行的优化就交给大家了：分块求平均。我眼睛花了。

本文参考