仿射变换大概是图像处理里最常见的一种操作了，旋转、缩放、平移、裁剪拼接，全都离不开它。我之前在做一个实时图像处理模块时，需要在一台普通电脑上同时处理多路视频帧，每帧都要做一次任意角度的旋转加缩放。CPU 上调用现成库的 warpAffine，单帧 1024x1024 灰度图就要吃掉十几毫秒，多路并跑根本扛不住。后来我把这段逻辑重写成 CUDA 核函数，直接在 GPU 上做反向映射和双线性插值，单帧耗时直接降到零点几毫秒量级，整条管线一下子松快了。这篇文章就把我从原理到实现的完整路径整理出来，包含矩阵求逆、像素中心约定、边界处理、常量内存和纹理内存加速这些细节，适合已经会写简单 CUDA kernel、但对图像几何变换实现细节还不够清楚的开发者，也适合正为图像预处理性能发愁的朋友参考。

1. 先拆明白仿射变换：数学公式、坐标方向与并行模型

1.1 仿射变换的数学本质：一张矩阵搞定旋转、缩放和平移

仿射变换本质上就是一个 2x3 矩阵作用在二维坐标上。假设源图像里的一个像素坐标是(srcX, srcY)，目标图像里对应点是(dstX, dstY)，那么它们的关系可以表示为：

dstX = a * srcX + b * srcY + tx dstY = c * srcX + d * srcY + ty

这里的a、b、c、d组成 2x2 的线性变换部分，负责旋转、缩放、错切；tx、ty是平移分量。把整个变换写成矩阵形式是：

| dstX | | a b tx | | srcX | | dstY | = | c d ty | * | srcY | | 1 | | 0 0 1 | | 1 |

旋转 30 度加缩放 0.8 再平移，对应的矩阵就是三个基础矩阵的乘积。以列向量约定为例，通常按"先缩放、再旋转、最后平移"的顺序组合：M = T * R * S。理解为最终矩阵先作用缩放，再作用旋转，最后加上平移。旋转矩阵长这样：

R = | cos(theta) -sin(theta) 0 | | sin(theta) cos(theta) 0 | | 0 0 1 |

缩放矩阵是：

S = | scaleX 0 0 | | 0 scaleY 0 | | 0 0 1 |

平移矩阵是：

T = | 1 0 tx | | 0 1 ty | | 0 0 1 |

三个矩阵乘起来，就得到最终的 2x3 参数。这里的顺序千万不能搞反。我见过不少同事把旋转和平移顺序写反，结果图像绕着错误中心旋转，要么转完位置跑偏，要么绕着原点转，整个画面直接飞出可视区。如果不想手动推公式，也可以先在 CPU 上用现成库把变换矩阵算好，再传给 CUDA，但务必确认坐标约定一致。

1.2 正向映射与反向映射：为什么 GPU 实现必须用反向

新手最容易犯的错误，是直接按正向映射写核函数：遍历源图像每个像素，通过矩阵算出它应该出现在目标图像的哪个位置，然后把灰度值写过去。这个思路本身没错，但在 GPU 上实现时会遇到两个致命问题。

第一个问题是目标图像会出现大量空洞。旋转后，某些目标像素位置上没有任何源像素落到这里，于是产生密集的黑点。第二次映射时，多个源像素可能映射到同一个目标位置，后者会覆盖前者，到底保留谁的结果？没法确定。解决正向映射空洞通常需要反向扩散或者多次迭代，这在 GPU 上成本很高。

第二个问题是写回时的竞争。多个线程可能写同一个目标位置，必须用原子操作或者加锁，性能立刻崩掉。

因此成熟的做法是反向映射。遍历目标图像的每个像素(dstX, dstY)，用逆矩阵算出它来自源图像的哪个位置(srcX, srcY)，再到源图像上采样灰度值。这样做有几个直接的好处：目标图像每个像素恰好被一个线程处理，完全没有写冲突；不会出现空洞，因为目标像素逐个填满；数据访问天然是只读源图像、只写目标图像，非常适合 GPU 并行模型。

反向映射的公式就是把正向矩阵求逆：

srcX = invA * dstX + invB * dstY + invTx srcY = invC * dstX + invD * dstY + invTy

其中invA、invB、invC、invD是 2x2 线性部分的逆，invTx、invTy是新的平移量。后面在主机端准备时会给完整求逆代码。

1.3 并行模型拆解：一个线程负责一个目标像素

CUDA 处理仿射变换的思路，和像素级图像处理完全一致：把目标图像划分成网格，每次 kernel launch 分配一个二维线程块网格，每个线程负责计算目标图像中一个像素的输出值。

线程与像素的对应关系通常这样定：线程threadIdx是块内编号，blockIdx是块编号，两者组合得到全局像素坐标：

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

如果x或y超出了目标图像尺寸，直接 return。这种写法叫线程块网格映射，是图像处理 kernel 的标准骨架。每个线程内部要干的事很简单：

1. 根据线程坐标得到目标像素坐标；
2. 用逆矩阵把目标坐标变换到源图像坐标；
3. 判断源坐标是否越界；
4. 在源图像上做双线性插值；
5. 把结果写回目标图像。

这个流程天然并行，因为第 i 行第 j 列的输出只依赖于源图像和变换矩阵，不依赖其他目标像素。也就是说，只要矩阵参数不变，任何线程都可以独立完成任务，GPU 上千个核心可以同时开工。

2. 核函数实现关键环节：坐标、插值与边界处理

2.1 主机端准备：矩阵合成顺序与求逆的计算

在写核函数之前，主机端要把变换矩阵合成并求逆。我习惯用一个结构体存仿射参数：

struct Warp2x3 { float a, b, tx; float c, d, ty; };

合成矩阵的时候按列向量约定写，顺序是M = T * R * S。假设我们想做"先缩放 1.2 倍、再旋转 45 度、最后平移 (10, 20)"，主机端代码是这样：

float scaleX = 1.2f, scaleY = 1.2f; float theta = 45.0f * 3.14159265f / 180.0f; float cosT = cosf(theta), sinT = sinf(theta); // 缩放矩阵 Warp2x3 S = { scaleX, 0.0f, 0.0f, 0.0f, scaleY, 0.0f }; // 旋转矩阵 Warp2x3 R = { cosT, -sinT, 0.0f, sinT, cosT, 0.0f }; // 平移矩阵 Warp2x3 T = { 1.0f, 0.0f, 10.0f, 0.0f, 1.0f, 20.0f };

两个 2x3 矩阵相乘，自己写个小函数。关键点是先缩放后旋转再平移，相当于T * R * S。如果你的参考实现是行向量约定，矩阵要转置，否则旋转方向会反，这一点很容易踩坑。

接下来求逆矩阵。2x2 线性部分的逆很好算，设：

det = a * d - b * c invA = d / det invB = -b / det invC = -c / det invD = a / det

平移部分的逆要乘上线性部分的逆：

Warp2x3 invert(const Warp2x3& m) { float det = m.a * m.d - m.b * m.c; Warp2x3 inv; inv.a = m.d / det; inv.b = -m.b / det; inv.c = -m.c / det; inv.d = m.a / det; inv.tx = -(inv.a * m.tx + inv.b * m.ty); inv.ty = -(inv.c * m.tx + inv.d * m.ty); return inv; }

求逆最好在 CPU 上用 double 算，算完再转成 float 传给 GPU。大图变换时 float 的累积误差会让边缘产生明显偏移，用 double 求逆能把这个误差压到最低。

2.2 像素中心约定：半像素偏移到底怎么来的

这是整个实现里最隐蔽、最容易出问题的地方。图像的每个像素其实是一个有面积的区域，我们平常说的像素坐标(i, j)到底指代像素的左上角，还是像素中心？大多数图像处理库采用像素中心约定，也就是整数坐标(i, j)对应第 i 行第 j 列像素的中心点。OpenCV 的 warpAffine 也是按这个约定处理坐标的，所以在 CPU 上调用看不出问题，但自己写核函数时就全体现在了坐标变换的偏移上。

按中心约定，目标图像像素(x, y)的中心坐标是(x + 0.5, y + 0.5)。用逆矩阵变换得到的源坐标是源图像里一个连续坐标，采样时也要用它代表源像素的中心。但如果你直接把变换后的坐标当作像素索引去取整，就会产生 0.5 像素的整体偏移，结果图像看起来像被推动过一样，细腻的边缘会变糊，缩放时尤其明显。

正确做法是全程使用中心坐标。代码里体现为：

float dstX = (float)x + 0.5f; float dstY = (float)y + 0.5f; float srcX = c_invM[0] * dstX + c_invM[1] * dstY + c_invM[2]; float srcY = c_invM[3] * dstX + c_invM[4] * dstY + c_invM[5]; // 采样前把源坐标从中心坐标转换回像素索引坐标 srcX -= 0.5f; srcY -= 0.5f;

最后那两行-0.5f作用是把源坐标从连续中心坐标转成以像素左上角为参考的索引坐标，这样floorf(srcX)得到的才是正确的左邻像素索引。整个过程相当于先加 0.5 再减 0.5，看似没有意义，但如果中间漏掉一步，最终结果一定偏半个像素。

2.3 双线性插值与边界策略：核函数主体代码

双线性插值的意思是，当变换后的源坐标落在四个像素之间时，按照距离权重把四个像素的值混合起来。公式不复杂，但实现时要注意浮点取整和边界判断的顺序。

假设源坐标为srcX、srcY，先取左邻和上邻的整数索引：

int x0 = (int)floorf(srcX); int y0 = (int)floorf(srcY); float fx = srcX - x0; float fy = srcY - y0;

fx、fy就是目标点在格子内的横向、纵向比例。用这四个像素和两个比例权重做两次水平插值、一次垂直插值即可。完整核函数长这样：

__constant__ float c_invM[6]; __global__ void affine_warp_linear( const unsigned char* __restrict__ src, unsigned char* __restrict__ dst, int srcW, int srcH, int dstW, int dstH) { int x = blockIdx.x * blockDim.x + threadIdx.x; int y = blockIdx.y * blockDim.y + threadIdx.y; if (x >= dstW || y >= dstH) return; float dstX = (float)x + 0.5f; float dstY = (float)y + 0.5f; float srcX = c_invM[0] * dstX + c_invM[1] * dstY + c_invM[2] - 0.5f; float srcY = c_invM[3] * dstX + c_invM[4] * dstY + c_invM[5] - 0.5f; int x0 = (int)floorf(srcX); int y0 = (int)floorf(srcY); float fx = srcX - x0; float fy = srcY - y0; if (x0 < 0 || y0 < 0 || x0 + 1 >= srcW || y0 + 1 >= srcH) { dst[y * dstW + x] = 0; return; } int x1 = x0 + 1; int y1 = y0 + 1; float v00 = (float)src[y0 * srcW + x0]; float v10 = (float)src[y0 * srcW + x1]; float v01 = (float)src[y1 * srcW + x0]; float v11 = (float)src[y1 * srcW + x1]; float top = (1.0f - fx) * v00 + fx * v10; float bot = (1.0f - fx) * v01 + fx * v11; float val = (1.0f - fy) * top + fy * bot; dst[y * dstW + x] = (unsigned char)(val + 0.5f); }

边界判断放在插值之前，是为了防止越界读。上面的写法是越界像素直接填 0，适合黑色背景场景。如果希望边界像素自然延伸，可以把越界处理改成钳制：

if (x0 < 0) x0 = 0; if (y0 < 0) y0 = 0; if (x0 + 1 >= srcW) x0 = srcW - 2; if (y0 + 1 >= srcH) y0 = srcH - 2;

钳制之后再用插值公式，效果是边缘像素向外重复，比填黑自然得多。两种方式都在生产环境里很常用，按需求取舍即可。

2.4 多通道扩展：RGB 图像怎么处理

上面的核函数只处理了单通道灰度图。实际项目中更多是三通道 RGB 或 RGBA。最简单的扩展方式是每个线程处理一个像素的三个通道，每个通道单独做双线性插值，然后写到目标图像的对应位置。可以把unsigned char换成uchar3，或者直接用float3在寄存器里做插值。

一个实用技巧是，三个通道可以复用一个线程内已经算好的srcX、srcY、fx、fy，也就是说坐标变换部分只算一次，插值部分重复三遍。这样写比把图像拆成三个单通道图、分别调用三次核函数快得多。CUDA 里uchar4是 128 位对齐访问，一次内存操作就能读入四个字节，比三次独立的 8 位读更友好。把源指针强转成const uchar4*后，一个线程处理一个像素的 BGRA，带宽利用率和访存效率都会更好。

3. 让核函数跑得更快：常量内存、纹理采样与访存调优

3.1 常量内存：6 个浮点参数的广播通道

仿射变换核函数里只用了 6 个浮点参数，这些参数在所有线程中完全相同，是常量内存的典型使用场景。用__constant__声明的数组会被放到 GPU 的常量内存空间，对所有线程广播访问时，速度比直接传全局内存指针或者每次都从参数里取值更快。

常量内存的声明和拷贝是这样的：

__constant__ float c_invM[6]; float h_invM[6] = { inv.a, inv.b, inv.tx, inv.c, inv.d, inv.ty }; cudaMemcpyToSymbol(c_invM, h_invM, sizeof(float) * 6);

需要注意，cudaMemcpyToSymbol只需要调用一次，只要变换参数不变，后续每次 launch 都不需要重新拷贝。如果每一帧的变换矩阵都不同，再考虑每次更新；但即便如此，6 个浮点的拷贝开销也几乎可以忽略。

这里有个小优化点：把 6 个参数放到__constant__里，还能减少 kernel 参数传递的寄存器占用，虽然效果微乎其微，但会让代码更清晰。老项目里也常见用一个__constant__ float c_mat[2][3]的二维数组，两者等价，按习惯来。

3.2 纹理内存：把插值和边界交给硬件

手写双线性插值已经有不错的性能，但还有更取巧的路径：使用纹理内存。GPU 的纹理采样单元天生支持双线性插值，而且硬件在纹理缓存上做了针对二维空间局部性的优化。旋转、缩放后的图像访问虽然整体是倾斜的，但相邻目标像素对应的源坐标仍然聚集在一起，纹理缓存能显著提高命中率。

用纹理对象实现仿射变换，核心优势有两个：一是省掉手写插值代码，核函数里一次tex2D调用就完成采样；二是边界处理可以交给纹理地址模式，使用cudaAddressModeClamp让越界坐标钳制到边缘像素，连边界分支判断都省了。

创建纹理对象的代码大致如下：

cudaTextureObject_t createTexObj(const float* d_src, int w, int h) { cudaResourceDesc resDesc = {}; resDesc.resType = cudaResourceTypePitch2D; resDesc.res.pitch2D.devPtr = (void*)d_src; resDesc.res.pitch2D.desc = cudaCreateChannelDesc<float>(); resDesc.res.pitch2D.width = w; resDesc.res.pitch2D.height = h; resDesc.res.pitch2D.pitchInBytes = w * sizeof(float); cudaTextureDesc texDesc = {}; texDesc.addressMode[0] = cudaAddressModeClamp; texDesc.addressMode[1] = cudaAddressModeClamp; texDesc.filterMode = cudaFilterModeLinear; texDesc.readMode = cudaReadModeElementType; texDesc.normalizedCoords = 0; cudaTextureObject_t tex = 0; cudaCreateTextureObject(&tex, &resDesc, &texDesc, nullptr); return tex; }

这里我把源图转成了float数组，是为了让线性采样模式更直接。实际项目中如果源图本来就是uchar，可以先用cudaChannelFormatKindUnsigned的通道格式创建纹理，但线性插值的精度和边界行为在不同显卡上有细微差别，我的经验是统一转成 float 会少很多麻烦。

核函数变得极短：

__global__ void affine_warp_tex( unsigned char* __restrict__ dst, int dstW, int dstH, cudaTextureObject_t tex) { int x = blockIdx.x * blockDim.x + threadIdx.x; int y = blockIdx.y * blockDim.y + threadIdx.y; if (x >= dstW || y >= dstH) return; float dstX = (float)x + 0.5f; float dstY = (float)y + 0.5f; float srcX = c_invM[0] * dstX + c_invM[1] * dstY + c_invM[2] - 0.5f; float srcY = c_invM[3] * dstX + c_invM[4] * dstY + c_invM[5] - 0.5f; float val = tex2D<float>(tex, srcX + 0.5f, srcY + 0.5f); dst[y * dstW + x] = (unsigned char)(val + 0.5f); }

tex2D里的+0.5f是因为纹理采样坐标 0.5 位置对应第一个像素的中心，而我上面的srcX已经减过 0.5 变成像素索引坐标，所以采样时要加回去。这个换算关系在不同资料里表述不一样，但只要始终盯着"像素中心在连续坐标里的位置"就不会错。

纹理版本的性能通常比手写插值版本快 30% 左右，代码量却少了一半。代价是纹理对象创建有一定开销，如果只是一次性变换，可能不划算；如果是在视频管线里反复做变换，纹理对象可以复用，收益巨大。

3.3 线程配置与访问优化：x 方向连续是第一原则

CUDA 里线程块尺寸对性能影响不小。我一般用dim3 block(16, 16)或dim3 block(32, 8)，网格按图像尺寸向上取整：

dim3 block(16, 16); dim3 grid((dstW + block.x - 1) / block.x, (dstH + block.y - 1) / block.y);

很多新手会把 block 设成(32, 32)甚至更大，结果性能不升反降。原因在于占用率和线程切换开销的平衡。16x16 的 block 占 256 个线程，属于比较稳妥的默认值，大多数显卡上都能跑到接近满占用。

访存方面，目标图像是按行存储的，y 方向相邻线程访问的行相距很远，x 方向相邻线程访问的地址却是连续的。所以应该让线程块在 x 方向跨度大于 y 方向跨度，比如 32x8 这个形状，就是让同一 warp 里的 32 个线程尽量落在同一行上，这样写回目标图像时一次内存事务能覆盖更多连续地址。

源图像读取的局部性没那么直观。旋转后的源坐标不再严格按行对齐，但相邻目标像素对应的源像素仍然聚集，所以纹理缓存或者 L2 缓存能接得住。如果用全局内存直接读源图，可以考虑用共享内存做局部缓存：把目标图像一个 tile 对应的源图像区域预先加载到共享内存。不过这个方法只适合旋转角度小、缩放比接近 1 的场景，角度一大，源区域就变成斜着的平行四边形，共享内存利用率下降，收益有限。我实际项目里更多依赖纹理缓存，代码简单且不容易出边界问题。

3.4 实测性能对比

我不建议背任何人的绝对数值，因为显卡型号、驱动、图像尺寸都会影响结果。但相对量级可以给你一个参考。假设 1024x1024 灰度图，参考设备是一块中端消费级显卡，数据已经在设备端，只算 kernel 耗时：

CPU 和 GPU 差距大约是 20 到 40 倍，纹理版本又在基础 CUDA 版本上快一点点。但要注意，如果数据要从主机内存拷贝到显存，一次 1MB 的拷贝大约 0.2 到 0.5 ms，可能比 kernel 本身还贵。这也是为什么我建议尽量把整条图像处理链路留在 GPU 上，而不是 CPU 和 GPU 来回倒腾。如果一定要拷贝，优先用页锁定内存或者统一虚拟内存来降低传输延迟。

4. 常见问题与排查技巧实录

4.1 图像错位半个像素：坐标约定的锅

症状是变换结果整体偏移 0.5 像素，放大图像时边缘发虚，旋转后细线出现台阶，但整体轮廓又像是对的。排查思路很简单：先用一个小尺寸图像（比如 4x4），在 CPU 上手动算几个点的坐标，跟 GPU 输出对比。只要你发现目标像素(0,0)应该对应源像素(0,0)的中心，但实际采样到的值来自(-0.5,-0.5)位置，问题就锁定在像素中心偏移上。

解决办法就是 2.2 节讲的那套流程：目标坐标加 0.5，变换后减 0.5。不少开源代码里能见到只用+0.5不用-0.5，或者反过来，那是因为他们在插值前用了不同的坐标基准。你只需要保证目标侧和源侧用同一套约定即可。

4.2 输出全黑或花屏：逆矩阵传错了

输出大部分区域是纯黑，只有个别像素有值，或者花屏一样的噪点，基本可以断定逆矩阵求错了。常见原因有三个：

第一个是直接把正向矩阵传给了核函数。正向往目标映射会计算出完全错误的源坐标，大部分坐标落在图像外，边界判断把它们全部填成黑，于是出现黑屏。排查方法是在主机端用同一个矩阵手算一个点的正向和反向映射，确认两者互逆。

第二个是矩阵元素拷贝顺序错乱。c_invM数组的索引顺序必须和核函数里的公式严格一致，我建议把数组声明和赋值注释写清楚，避免索引错位。

第三个是矩阵组合顺序不对。先旋转再平移和先平移再旋转会产生完全不同的结果，如果你的图像旋转后位置绕着左上角转，而不是绕着图像中心转，多半就是你忘了先把坐标系移到中心。处理图像旋转时，通常要先把图像中心平移到原点，旋转完再平移回去，等价于把平移分量修改为txNew = tx + cx - cosT * cx + sinT * cy这样的形式。这个绕中心的修正，是仿射变换里最经典的一个坑。

4.3 性能上不去：边界分支和访存问题

如果你发现核函数跑得很慢，先看占用率。边界分支判断本身开销不大，但写得不好的话，每个线程都要跑一堆if，指令数翻倍。用纹理版本加 Clamp 地址模式，可以把边界分支完全消掉，这是立竿见影的优化。

另一个性能问题是源图像全局内存访问不连续。旋转后，线程读取的源坐标呈倾斜排列，同一 warp 里的线程访问的地址跨度很大，导致每次内存事务都包含大量无用数据。这时候纹理缓存的收益就体现出来了。如果在普通全局内存上硬读，可以试着把 block 尺寸从 16x16 改成 32x8，让同一个 warp 的那 32 个线程尽量靠近，减少行跨度。

最后，别忘了检查有没有因为索引计算越界导致的非法访问，CUDA 的非法访问会把整个 kernel 崩掉，日志里能看到out of bounds相关报错。先用compute-sanitizer跑一遍小图，能快速定位越界位置。

4.4 常见问题速查表

这个速查表基本覆盖了我在项目里踩过的大部分坑。如果你写的变换矩阵还带了错切或者不规则缩放，建议额外在主机端写一个单点验证函数，把几个关键像素的变换结果打印出来，和 CPU 参考输出对比，排查起来会快很多。

说实话，第一次跑通 CUDA 仿射变换的时候，我盯着输出图像愣了半天。同样一个变换，CPU 上循环上百万次，GPU 上就是一次 kernel launch 的事。后来做性能优化时，我又花了不少时间在纹理内存和坐标约定上，才发现很多看似玄学的图像错位，根源就是半像素偏移。这里分享几个经验：主机端永远留一个 CPU 参考实现，方便单点验证；边界处理决定上线后会不会被截图吐槽，Clamp 通常比填黑更安全；纹理内存在大角度旋转加缩放的场景下收益最明显，值得优先尝试。如果后续有时间，我再聊聊用共享内存把局部性吃透的做法，以及多路视频并行时的流管理细节。