第12章：CUDA 图像处理：内核与内存模型

第 11 章装好了环境，这一章让 GPU 真正开始处理图像。从第一性原理出发：一条 for 循环如何变成几十万个线程、线程如何编组成 grid/block/thread、数据如何住在 global/shared/寄存器三级内存里——最后亲手编译运行两个真正的图像内核：像素翻转与灰度转换

12.1 第一性原理：从循环到内核

上一章我们把 CUDA 环境立了起来：驱动在跑、nvcc 就位、deviceQuery 打出 Result = PASS。工具备好，这一章开始真正的 CUDA 图像处理。先亮出方法论——费曼的第一性原理：不问"CUDA 的 API 怎么背"，先问"图像处理最底层的运算是什么"。答案朴素得很：对每个像素做同一件独立的事——取反、变灰、找边缘。第 11 章 11.1 节已经算过这笔账：一张 1920×1080 的图约 207 万像素，逐像素独立运算，天然适合并行。CPU 版是一条 for 循环串行扫 200 万次，GPU 版是把循环体交给 200 万个线程同时执行。第一性原理的推论水到渠成：把"循环体"变成"内核函数体"，把"循环变量"变成"线程编号"，一条 for 循环就翻译成了 CUDA 内核。

官方编程指南（CUDA C++ Programming Guide）对内核（kernel）的定义很朴素：在 GPU 上执行的函数就叫内核，启动内核（launch）相当于同时启动成百上千万个线程并行执行同一段代码。语法上内核有两条硬规矩：① 用 __global__ 修饰——告诉编译器这个函数在设备端执行、由主机端启动；② 返回类型必须是 void——官方明说，内核无法直接返回值，计算结果只能写进全局内存（显存），再由主机用 CUDA API 拷回。还要先建立 host/device 这对概念：host 指 CPU 及其内存，device 指 GPU 及其显存（第 11 章提过"异构计算"，这里给出官方定义）；程序永远从 CPU 启动，主机代码负责拷贝数据、启动内核、等待完成。下面把第 11 章 11.1 节那个"像素取反"的伪代码变成第一个真内核。

// invert_kernel.cu —— 第一个真内核：灰度图像素取反
__global__ void invert(unsigned char* data, int n) {
    // 一维索引：本线程在整个网格中的全局编号（12.2 节详解）
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {               // 边界检查：网格可能比数据大
        data[i] = 255 - data[i]; // 取反：亮度 255 变 0，0 变 255
    }
}

// kernel_min.cu —— 最小可运行 CUDA 程序：一维数组取反
#include <stdio.h>
#include <cuda_runtime.h>

__global__ void invert(unsigned char* data, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) data[i] = 255 - data[i];
}

int main() {
    const int N = 1024;                       // 一维"图像"：1024 个像素
    unsigned char* h = (unsigned char*)malloc(N);
    for (int i = 0; i < N; ++i) h[i] = i % 256; // 造一条渐变数据

    unsigned char* d = NULL;
    cudaMalloc(&d, N);                         // 1. 在显存分配
    cudaMemcpy(d, h, N, cudaMemcpyHostToDevice);   // 2. 拷入显存
    invert<<<1, 256>>>(d, N);               // 3. 启动内核：1 个块 256 线程
    cudaMemcpy(h, d, N, cudaMemcpyDeviceToHost);   // 4. 结果拷回内存
    cudaFree(d);                                   // 5. 释放显存

    printf("h[0]=%d h[255]=%d (期望 255, 0)\n", h[0], h[255]);
    free(h);
    return 0;
}

对照第 11 章 11.1 节的两个代码块：CPU 版把"取反"写进两层 for 循环，GPU 版把它写进 invert 函数体；CPU 版的循环边界是 y < H、x < W，GPU 版的"循环边界"变成了 blockIdx * blockDim + threadIdx 与一个 if 守卫。核心思想一句话：循环是串行的枚举，内核是并行的枚举。剩下的问题全是工程细节：线程怎么编组（12.2）、像素坐标怎么算（12.3）、数据放在哪（12.4）、怎么编译运行（12.6）。

12.2 线程层次：grid、block 与 thread

内核一启动就是几十万上百万个线程，它们不是散兵游勇，而是有严格的编组——官方叫线程层次（thread hierarchy）：线程（thread）组成线程块（thread block），线程块组成网格（grid），网格与块都可以是一维、二维或三维。为什么中间要夹一层"块"？这是从硬件事实推出来的（第一性原理）：GPU 由若干流式多处理器（SM）组成，一个线程块的所有线程被调度到同一个 SM 上执行，块是调度的基本单位；而 GPU 上的 SM 只有几十上百个，网格里的块可能有几十万个，所以官方明确要求：块之间没有执行顺序保证、不得互相依赖——这样任意块都可以按任意顺序塞进任意 SM。推论很重要：块内可以协作（共享内存、同步），块间必须各自独立，这决定了内核的写法。

每个线程靠四个内置变量认清自己是谁：threadIdx（本线程在块内的编号）、blockIdx（本块在网格中的编号）、blockDim（块的大小，每维多少线程）、gridDim（网格的大小，每维多少块），各自带 .x/.y/.z 三个分量。一维配置下，线程的全局唯一编号 = blockIdx.x * blockDim.x + threadIdx.x。还有一个硬件细节要记住：线程按 32 个一组组成 warp，同一 warp 的线程按 SIMT 模式齐步执行指令（官方术语），所以块大小最好取 32 的倍数（比如 256 = 8 个 warp），否则最后一个 warp 会空转；官方还给出硬上限：每块最多 1024 个线程（设备属性 maxThreadsPerBlock）。下面用两个小程序把层次"看"出来。

// whoami.cu —— 每个线程报出自己的"住址"：块号、块内号、全局号
#include <stdio.h>
#include <cuda_runtime.h>

__global__ void whoami() {
    int block = blockIdx.x;     // 本线程所在块：0 .. gridDim.x-1
    int thread = threadIdx.x;   // 块内线程号：0 .. blockDim.x-1
    int global = block * blockDim.x + thread; // 全局唯一编号
    printf("block %d, thread %d -> global %d\n", block, thread, global);
}

int main() {
    whoami<<<2, 4>>>();      // 2 个块、每块 4 个线程，共 8 个线程
    cudaDeviceSynchronize(); // 内核异步执行，等它跑完再退出（12.6 节详讲）
    return 0;
}

// vec_invert.cu —— 一维大数组：网格大小按数据量计算
__global__ void invert(unsigned char* data, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) data[i] = 255 - data[i];
}

int main() {
    const int N = 1 << 20;        // 1M 个像素（约 100 万）
    int blockSize = 256;          // 每块 256 线程 = 8 个 warp（32 的倍数）
    int gridSize = (N + blockSize - 1) / blockSize; // 上取整：4096 块
    // ... 分配与拷贝略（12.4 节给完整五步流程）
    invert<<<gridSize, blockSize>>>(d, N);
    // ... 回拷与验证略
    return 0;
}

注意 vec_invert 里网格大小的算法：(N + blockSize - 1) / blockSize 是"上取整"——N 除以 blockSize 有余数时，也得凑够一个块，多出来的线程靠 if 守卫跳过。块数远超 SM 数没关系：SM 跑完一个块就接下一个，官方保证"任意大的网格都能在任意小的 GPU 上跑"。不过对图像来说，一维编号还不够直观——像素天生是二维的，下一节把网格也变成二维。

12.3 二维索引：把像素坐标算出来

图像是二维数据：宽 W、高 H，在内存里按行优先排成一维数组（第 9 章讲 Mat 时就是这个布局，这里复用）。要把"线程"和"像素"对上号，最自然的做法是让网格和块都取二维：块里 16×16 个线程，网格里若干块铺满整张图。每个线程的像素坐标由两行公式算出：x = blockIdx.x * blockDim.x + threadIdx.x、y = blockIdx.y * blockDim.y + threadIdx.y，再折算回一维下标 idx = y * width + x——这正是第 11 章 11.1 节伪代码里那两个公式的正式版。官方编程指南明确说，网格与块支持多维，目的就是方便把线程映射到二维数据上。

二维映射立刻暴露一个新手必踩的坑：图像尺寸往往不是块大小的整数倍。640×480 的图配 16×16 的块，网格上取整后是 40×30 块，多出来的线程对应"不存在的像素"——如果放任不管，就会越界读写显存，轻则数据损坏、重则崩溃（官方术语叫未定义行为）。标准解法是边界守卫：if (x < width && y < height)，越界的线程直接跳过。启动配置里网格同样用上取整公式：dim3 grid((W + 15) / 16, (H + 15) / 16)。下面先写带边界检查的二维取反内核，再写一个"画棋盘"的内核——把 (x + y) 的奇偶映射成黑白，跑完直接得到一张肉眼可见的 PGM 图片，索引算得对不对一目了然。

// invert2d.cu —— 二维网格/块：一个线程负责一个像素（含边界守卫）
__global__ void invert(unsigned char* img, int width, int height) {
    int x = blockIdx.x * blockDim.x + threadIdx.x; // 像素列
    int y = blockIdx.y * blockDim.y + threadIdx.y; // 像素行
    if (x < width && y < height) {     // 越界守卫
        int idx = y * width + x;     // 行优先：二维坐标折算成一维下标
        img[idx] = 255 - img[idx];
    }
}

int main() {
    const int W = 640, H = 480;
    // ... 造图、cudaMalloc、cudaMemcpy 略（12.5 节给完整程序）
    dim3 block(16, 16);                    // 每块 16x16 = 256 线程（8 个 warp）
    dim3 grid((W + 15) / 16, (H + 15) / 16); // 上取整：40 x 30 = 1200 块
    invert<<<grid, block>>>(d_img, W, H);
    // ... 回拷、验证、释放略
    return 0;
}

// checker.cu —— 用 (x + y) 的奇偶生成棋盘并写成 PGM，肉眼验证二维索引
#include <stdio.h>
#include <cuda_runtime.h>

__global__ void checker(unsigned char* img, int width, int height) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if (x < width && y < height) {
        int idx = y * width + x;
        img[idx] = ((x + y) % 2 == 0) ? 255 : 0; // 奇偶交替：黑一格白一格
    }
}

int main() {
    const int W = 256, H = 256;
    unsigned char* h = (unsigned char*)malloc(W * H);
    unsigned char* d = NULL;
    cudaMalloc(&d, W * H);

    dim3 block(16, 16);
    dim3 grid((W + 15) / 16, (H + 15) / 16); // 16 x 16 = 256 块
    checker<<<grid, block>>>(d, W, H);
    cudaMemcpy(h, d, W * H, cudaMemcpyDeviceToHost);

    FILE* f = fopen("checker.pgm", "wb"); // PGM P5：文本头 + 二进制像素
    fprintf(f, "P5\n%d %d\n255\n", W, H);
    fwrite(h, 1, W * H, f);
    fclose(f);
    cudaFree(d);
    free(h);
    return 0;
}

跑完 checker.cu 用图片查看器打开 checker.pgm，应该是一张 256×256 的黑白棋盘。如果棋盘错位、花屏，多半是下标折算错了——这就是"可见的调试"：用图案验证索引，比打印一万个数字直观得多。二维索引掌握后，任意"逐像素独立"的运算（取反、调亮度、加噪声）都能照这个模板写。但还有一个大问题没解决：像素数据在内存里，GPU 怎么拿到？下一节讲内存层次与数据搬运。

12.4 内存层次：global、shared 与寄存器

数据放在哪、多快能取到，是 GPU 编程的命门。官方内存模型把 GPU 的内存分成几层——第一性原理的推论：内存离计算单元越近就越快、也越小。从快到慢排：寄存器（register）在 SM 片上，线程私有，编译器自动分配，放线程的局部变量，最快；共享内存（shared）也在 SM 片上，块内所有线程可见，比全局快一个数量级，官方说它与 L1 缓存共用同一块物理资源；全局内存（global）是 GPU 板载 DRAM，整个网格所有线程都能访问，容量最大但延迟最高，官方原话：全局内存的分配持续到程序结束或调用 cudaFree。官方内存类型表还列出另外两种：常量内存（constant），网格可见的只读数据，典型 64KB，适合所有人共用的参数；局部内存（local），逻辑上线程私有、物理上却在显存——寄存器不够用时编译器会把变量"溢出"到这里，访问它跟访问全局一样慢。

内存类型	作用域	生命周期	物理位置	速度（相对）
寄存器 register	线程私有	内核执行期	SM 片上	最快
共享 shared	块内共享	内核执行期	SM 片上（与 L1 同资源）	快
全局 global	网格可见	程序结束 / cudaFree	板载 DRAM	慢
常量 constant	网格可见（只读）	程序结束	板载 DRAM	慢（带缓存）
局部 local	线程私有	内核执行期	板载 DRAM（逻辑私有）	慢（寄存器溢出的后备）

host 与 device 内存物理分离，数据进出 GPU 全靠三个运行时 API（官方文档逐一核实）：cudaMalloc(&devPtr, size) 在显存分配 size 字节线性内存，把指针写进 devPtr（官方原话：分配的内存不清零）；cudaMemcpy(dst, src, count, kind) 拷贝 count 字节，第四个参数 kind 指定方向——cudaMemcpyHostToDevice（内存到显存）、cudaMemcpyDeviceToHost（显存到内存），另有 HostToHost 与 DeviceToDevice；cudaFree(devPtr) 释放，与 cudaMalloc 一一对应。标准流程五步：分配 → 拷入 → 启动内核 → 拷回 → 释放，下面的 mem_flow 就是这个骨架。共享内存用 __shared__ 声明（数组大小编译期定死），块内线程共享，配合 __syncthreads() 做块内同步——官方语义：阻塞块内所有线程，直到大家都执行到这一行，防止"我还没写你就来读"的竞态。

// mem_flow.cu —— 显存生命周期：分配 -> 拷贝 -> 内核 -> 回拷 -> 释放
int main() {
    const int BYTES = 640 * 480;        // 一张 640x480 灰度图
    unsigned char* h_in = (unsigned char*)malloc(BYTES);
    unsigned char* h_out = (unsigned char*)malloc(BYTES);
    // ... 读图/造图省略（12.5 节给完整程序）

    unsigned char* d_in = NULL;
    unsigned char* d_out = NULL;
    cudaMalloc(&d_in, BYTES);            // 分配显存（不清零，官方原话）
    cudaMalloc(&d_out, BYTES);
    cudaMemcpy(d_in, h_in, BYTES, cudaMemcpyHostToDevice);   // 内存 -> 显存
    kernel<<<grid, block>>>(d_in, d_out, W, H);  // 内核读写显存
    cudaMemcpy(h_out, d_out, BYTES, cudaMemcpyDeviceToHost); // 显存 -> 内存
    cudaFree(d_in);                      // 释放，与 cudaMalloc 一一对应
    cudaFree(d_out);
    // ... 验证 h_out 并收尾
}

// shared_demo.cu —— __shared__：块内共享内存（块内可见、生命周期同内核）
__global__ void tile_invert(unsigned char* img, int width, int height) {
    __shared__ unsigned char tile[16][16]; // 静态声明，编译期定大小

    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    int idx = y * width + x;

    // 1. 全体线程协作：把各自的像素读进共享内存
    if (x < width && y < height) {
        tile[threadIdx.y][threadIdx.x] = img[idx];
    }
    __syncthreads(); // 2. 同步：等块内所有线程写完 tile 再继续（官方语义）

    // 3. 从共享内存读（快），取反后写回全局内存
    if (x < width && y < height) {
        unsigned char v = tile[threadIdx.y][threadIdx.x]; // v 放寄存器
        img[idx] = 255 - v;
    }
}

老实说，这个 tile_invert 用共享内存并没有提速——每个像素只被读一次，直接读写全局就够了。它的价值是展示语法与心智模型：共享内存是"块内协作的工作台"，寄存器是"你手上的便签"，全局内存是"大家共用的仓库"。什么时候真正需要工作台？当多个线程反复用同一批数据、或块内要先算局部结果再合并时——这正是第 15 章性能优化方法论的系统主题，这里先埋个钩子。

12.5 图像处理内核实战：像素翻转与灰度转换

理论齐了，动手。第一个内核：灰度图像素翻转（取反），公式 out = 255 - in，对每个像素独立成立，负片效果。完整程序走一遍五步流程：主机端造一张 640×480 的渐变测试图 → cudaMalloc 分配显存 → cudaMemcpy 拷入 → 二维网格启动内核 → cudaDeviceSynchronize 等内核跑完（内核启动是异步的，官方明说 launch 立即返回，回拷前必须等它完成）→ 拷回验证 → 释放。运行后打印两个像素：h[0] = 255（原值 0 取反）、h[100] = 155（原值 100 取反），即成功。

// invert.cu —— 完整程序：灰度图取反（合成渐变测试图 + 打印验证）
#include <stdio.h>
#include <cuda_runtime.h>

__global__ void invert(unsigned char* img, int width, int height) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if (x < width && y < height) {
        img[y * width + x] = 255 - img[y * width + x];
    }
}

int main() {
    const int W = 640, H = 480;
    unsigned char* h = (unsigned char*)malloc(W * H);
    for (int i = 0; i < W * H; ++i) h[i] = i % 256; // 横向渐变

    unsigned char* d = NULL;
    cudaMalloc(&d, W * H);
    cudaMemcpy(d, h, W * H, cudaMemcpyHostToDevice);

    dim3 block(16, 16);
    dim3 grid((W + 15) / 16, (H + 15) / 16);
    invert<<<grid, block>>>(d, W, H);
    cudaDeviceSynchronize();            // 等内核跑完（启动是异步的）
    cudaMemcpy(h, d, W * H, cudaMemcpyDeviceToHost);

    printf("h[0]=%d h[100]=%d (期望 255, 155)\n", h[0], h[100]);
    cudaFree(d);
    free(h);
    return 0;
}

第二个内核：彩色转灰度。RGB 图每个像素占 3 个字节（本节按 R,G,B 连续排的自定义布局；第 9 章 OpenCV 的 Mat 是 BGR 顺序——本质一样，只是通道顺序不同，接 OpenCV 时要注意，这是它的老传统），灰度图每像素 1 字节。转换公式用 BT.601 亮度公式：Y = 0.299R + 0.587G + 0.114B——OpenCV 的 cvtColor 用的就是这套系数，呼应第 9、10 章。内核里每个线程负责一个像素：用二维索引算出像素号 px 与通道偏移 off = px * 3，读出三个通道，加权求和写进灰度数组。注意：这个内核同样逐像素独立、不需要共享内存——什么时候该上共享内存、块内协作怎么组织，留给第 15 章的性能方法论。

// rgb2gray.cu —— 彩色转灰度：每个线程处理一个像素的 3 个字节
// 布局：按像素连续存 RGB（R,G,B,R,G,B,...），灰度图每像素 1 字节
__global__ void rgb2gray(const unsigned char* rgb, unsigned char* gray,
                          int width, int height) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if (x < width && y < height) {
        int px = y * width + x;    // 第 px 个像素
        int off = px * 3;         // 三个通道的起始字节
        unsigned char r = rgb[off + 0];
        unsigned char g = rgb[off + 1];
        unsigned char b = rgb[off + 2];
        // BT.601 亮度公式（OpenCV cvtColor 同款系数）：
        gray[px] = (unsigned char)(0.299f * r + 0.587f * g + 0.114f * b);
    }
}

// 启动配置与 invert.cu 完全一样：
//   dim3 block(16, 16);  dim3 grid((W + 15) / 16, (H + 15) / 16);
//   rgb2gray<<<grid, block>>>(d_rgb, d_gray, W, H);

两个内核一灰一彩，覆盖了图像处理最基础的两种"逐像素映射"：单字节变换（取反）与多字节聚合（三通道加权）。它们都用同样的骨架：二维索引 + 边界守卫 + 一个公式。记住这个骨架，本系列的绝大多数图像算子（亮度、对比度、阈值、噪声）都能往里套。接下来的问题是：写好了、怎么变成能跑的程序？

12.6 编译运行：nvcc 与启动配置

编译靠 nvcc——第 11 章装好的 Toolkit 编译器，一条命令出可执行文件：nvcc -o invert invert.cu。它内部是两段式的：先调宿主 gcc 编译 CPU 部分，再用自己的编译器把内核编成 GPU 代码，最后链接成一个二进制——所以你只需要记住 nvcc 一个命令。运行 ./invert 直接看输出。如果报 nvcc: command not found，回到第 11 章检查 PATH（工具早已备好，别让环境问题挡住实验）。常用选项：-O2 开优化、-arch=sm_XX 指定目标架构（如 -arch=sm_86 对应 Ampere 代显卡；不指定也能跑，只是没吃满新架构特性）。

# 编译：nvcc 一条命令完成 host 部分（调 gcc）与 device 部分的编译链接
nvcc -o invert invert.cu

# 运行：内核异步执行，但程序里的 cudaMemcpy 回拷会隐式等内核完成
./invert
# 输出示例：
#   h[0]=255 h[100]=155 (期望 255, 155)

# 常用选项：-arch 指定目标架构，-O2 优化（细节见 nvcc 官方文档）
nvcc -arch=sm_86 -O2 -o invert invert.cu

启动配置 <<<grid, block>>> 是内核调用里最容易写错的部分，拆开看：第一参数是网格维度，第二参数是块维度，类型是内置向量类型 dim3——只写一个数时等价于一维（dim3 block(16, 16) 是二维块，没写的那维默认 1）。两条硬约束再强调一遍：块内线程数不超过 1024，且最好是 32 的倍数。还有个新手最迷惑的点：内核启动是异步的，launch 那一行本身不报执行错误——启动配置非法（比如块太大）也只会在之后暴露。所以官方推荐用 cudaGetLastError() 查启动是否合法、用 cudaDeviceSynchronize() 等内核结束并捕获执行期错误。把它们包进一个 CHECK 函数，每步都查——这是 CUDA 工程的基本功，也是后面排错的前提（编译报错、链接失败怎么解，第 13 章系统讲；运行期崩溃怎么调，第 14 章调试实战）。

// check_cuda.c —— 把错误检查包成函数：任何 CUDA 调用失败立即报错退出
void check(cudaError_t err, const char* file, int line) {
    if (err != cudaSuccess) {
        fprintf(stderr, "CUDA error %s at %s:%d\n",
                cudaGetErrorString(err), file, line);
        exit(1);
    }
}

#define CHECK(call) check((call), __FILE__, __LINE__)

// 用法：
//   CHECK(cudaMalloc(&d, N));           // 分配失败立刻报错
//   invert<<<grid, block>>>(d, W, H); // 内核启动本身不报错
//   CHECK(cudaGetLastError());          // 官方推荐：查启动配置是否合法
//   CHECK(cudaDeviceSynchronize());     // 等内核结束，查执行期错误

这一章我们从"一个像素的独立运算"出发（12.1），推导出线程的 grid/block/thread 编组（12.2）、二维像素索引（12.3）、global/shared/寄存器三级内存与数据搬运（12.4），并亲手编译运行了像素翻转与灰度转换两个图像内核（12.5、12.6）。环境装好了、内核跑通了，接下来的问题开始浮出水面：编译报错找不到符号怎么办？——第 13 章讲依赖与链接排错；程序崩溃怎么定位？——第 14 章调试实战；想让内核更快？共享内存、缓存友好与多线程的性能方法论，第 15 章系统展开。第一性原理的收益就在这里：你不再背 API，而是从"像素与独立运算"出发，自己推导出 CUDA 的每一块拼图。