MATLAB图像旋转：三种插值算法原理与手写实现

做数字图像处理实验的时候图像旋转这个题目几乎是绕不开的。很多人在MATLAB里直接一句imrotate就把结果拿到了但等到“不用内置函数手写实现”的要求一出来才发现里面名堂不少。最近也有好几个人问起“MATLAB实现图像旋转的三种算法”这个课设题我就把常用的三种插值旋转算法——最近邻、双线性、双三次——连原理带完整代码一起梳理一遍既适合交实验报告也适合真正想把旋转这块原理吃透的人。先把话说明白图像旋转里说的“三种算法”常见的理解有两种。一种是从坐标变换角度分有正向映射、逆向映射、三步错切分解另一种是从像素取值策略分也就是最近邻插值、双线性插值、双三次插值。实际做课设和项目时绝大多数语境下“三种算法”指的是后者。这篇博文以三种插值算法为主线但坐标映射那块是整个实现的地基我会把正向映射和逆向映射的坑也一并讲清楚。1. 从原理讲起图像旋转的本质与算法分类1.1 二维旋转的数学基础图像旋转本质上是对图像上每一个像素点做一次二维平面内的坐标变换。假设图像是灰度图像素在平面上的位置用坐标(x, y)表示绕原点逆时针旋转角度θ变换公式是x x·cosθ - y·sinθ y x·sinθ y·cosθ写成矩阵形式就是[x] [cosθ -sinθ] [x] [y] [sinθ cosθ] [y]这个矩阵就是二维旋转矩阵。但图像里有个特殊情况需要注意图像坐标系的y轴是向下的也就是像素的行号越往下越大这和数学里常用的直角坐标系y轴向上正好相反。所以在图像坐标系里直接套这个旋转矩阵角度方向会跟视觉直觉相反。MATLAB自带的imrotate以逆时针为正方向我也是按这个约定来写的。这样做的好处很明显方便和imrotate的结果做对照。你在写自己的旋转函数时最好也统一成这个方向不然后续对比调试容易越搞越乱。1.2 正向映射与逆向映射为什么不直接算确定用旋转矩阵之后第一个想到的做法通常是“正向映射”遍历原图的每个像素点用上面的公式算它旋转后的新位置把灰度值搬过去。问题来了。原图是离散网格旋转后坐标几乎不可能刚好落在整数位置上。比如原图像素(100, 80)绕中心旋转30度算出来的新坐标可能是(143.62, 51.37)那这个灰度值应该填到新图像的哪个像素直接四舍五入到(144, 51)原本多个原图像素旋转后可能同时映射到同一个目标像素产生重叠而另一些目标位置没有人映射到就成了空洞图像上会出现一条条规律的黑色纹路。所以工程上几乎不用正向映射直接做旋转都改用逆向映射遍历目标图像的每一个整数像素坐标用旋转矩阵的逆矩阵把它映射回原图坐标系。这样每个目标像素都能且只能拿到一个原图坐标从根源上避免了重叠和空洞问题。逆变换公式就是正向公式的角度取负x x·cosθ y·sinθ y -x·sinθ y·cosθ1.3 三种插值算法到底在插什么逆向映射得到的原图坐标一般也不是整数。比如目标像素(120, 80)映射回原图得到(98.36, 65.78)那这个位置的灰度值是多少原图上根本没有这个像素点只有(98, 65)、(99, 65)、(98, 66)、(99, 66)四个整数位置有像素值。“根据周围已知像素猜出中间位置像素值”的过程就叫插值。三种算法的差异核心就是“用多少个邻居、怎么加权”的区别最近邻插值找离映射坐标最近的整数像素直接抄它的值。邻居数量1个。双线性插值用周围4个像素做加权平均权重由距离决定。邻居数量4个。双三次插值用周围16个像素做加权平均权重由三次卷积核函数算出。邻居数量16个。一句话概括用的邻居越多、加权方式越讲究结果越平滑细腻但计算量也越大。这三种算法在实际工程里各有各的应用场景没有绝对的好坏。2. MATLAB手写三种插值旋转算法2.1 整体框架画布扩张与坐标逆映射动手写代码之前先把骨架搭好。整体流程分四步读入灰度图转成double类型方便做浮点运算。根据旋转角度计算旋转后的画布尺寸目标图的宽高和原图不一样。构造目标图像的坐标网格逐点做逆向映射得到每个目标像素对应的原图坐标。根据选定的插值方式在映射坐标处取灰度值。画布尺寸的计算有个简单通用的方法把原图四个角点绕图像中心旋转得到旋转后四个角点的坐标然后取x方向的最大最小值之差作为新宽度y方向的最大最小值之差作为新高度。这样连公式都不用背任何角度都不会出错。角度为θ时旋转后的包围盒尺寸也可以直接算rot_w ceil(|w·cosθ| |h·sinθ|) rot_h ceil(|w·sinθ| |h·cosθ|)这里w和h是原图宽和高。注意到没有新宽是由原宽乘以cos的绝对值加上原高乘以sin的绝对值构成的。90度时会退化成rot_w h、rot_h w和常识吻合。主函数写出来是这样function rot_img my_rotate(img, angle_deg, method) % 灰度图旋转支持三种插值方式 % img: 输入灰度图uint8 或 double 均可 % angle_deg: 旋转角度逆时针为正 % method: nearest / bilinear / bicubic % 统一转成 double 灰度图数值范围 [0, 1] if ~isa(img, double) img double(img) / 255; end [h, w] size(img); theta angle_deg * pi / 180; ct cos(theta); st sin(theta); % 旋转后画布尺寸 rot_h ceil(abs(w * st) abs(h * ct)); rot_w ceil(abs(w * ct) abs(h * st)); % 目标图像坐标系 [X, Y] meshgrid(1:rot_w, 1:rot_h); % 图像中心的坐标注意像素中心是 (size1)/2不是 size/2 cx1 (w 1) / 2; cy1 (h 1) / 2; cx2 (rot_w 1) / 2; cy2 (rot_h 1) / 2; % 逆向映射目标坐标 - 原图坐标 Xq ct .* (X - cx2) - st .* (Y - cy2) cx1; Yq st .* (X - cx2) ct .* (Y - cy2) cy1; % 三种插值 switch lower(method) case nearest rot_img nearest_interp(img, Xq, Yq); case bilinear rot_img bilinear_interp(img, Xq, Yq); case bicubic rot_img bicubic_interp(img, Xq, Yq); otherwise error(method 必须是 nearest、bilinear 或 bicubic); end % 截断到 [0,1] 并转回 uint8 rot_img max(0, min(1, rot_img)); rot_img uint8(round(rot_img * 255)); end这段代码里有个细节我得专门提一下图像中心的坐标用的是(size1)/2不是size/2。比如一张宽度为256的图像素横坐标从1到256中心和两边等距的位置是(2561)/2128.5不是128。很多人在这一步用size/2看起来差不多实际旋转后会整体偏移半个像素后面我专门写一节来排查这个问题。2.2 最近邻插值最朴素也最快最近邻插值很好理解映射回去的原图坐标是小数那就四舍五入取最近的整数像素直接赋值。function out nearest_interp(img, Xq, Yq) % 最近邻插值 [h, w] size(img); xr round(Xq); yr round(Yq); % 判断坐标是否在有效范围内 valid xr 1 xr w yr 1 yr h; % 把越界坐标 clamp 到边界内方便统一取数 xr_c min(max(xr, 1), w); yr_c min(max(yr, 1), h); out zeros(size(Xq)); idx sub2ind([h, w], yr_c, xr_c); out(valid) img(idx); end这里的sub2ind用来把行列坐标转成线性索引。注意它的参数顺序sub2ind([h, w], row, col)先是行号再是列号。我在代码里yr_c是行方向xr_c是列方向所以顺序不能反否则取出来的像素位置就串了。没有通过valid判断的位置保持默认值0也就是黑色背景。旋转后四个角上的三角形区域没有原图像素覆盖显示为黑色这是常规做法。最近邻的视觉效果可以用四个字形容颗粒感强。旋转后斜边会出现明显的锯齿像楼梯一样一格一格放大看非常明显。但它的优点是速度极快因为只做了一次四舍五入和一次查表没有乘加浮点运算。如果只是做旋转后的尺寸预览、缩略图生成或者对速度要求极高而质量要求不高的场景最近邻是可以考虑的。2.3 双线性插值平滑视觉的主力双线性插值的思想是用距离映射坐标最近的4个像素按距离远近做加权平均。位置左上、右上、左下、右下四个网格点的灰度值分别记为I00、I01、I10、I11映射点到左上角的水平和垂直距离分别为dx、dy则输出值out (1-dx)·(1-dy)·I00 dx·(1-dy)·I01 (1-dx)·dy·I10 dx·dy·I11这个公式的权重很有意思离哪个像素越近那个像素的权重越大。当dx和dy都为0时结果就是左上角像素本身完美还原当dx和dy都为1时结果就是右下角像素本身。function out bilinear_interp(img, Xq, Yq) % 双线性插值 [h, w] size(img); x0 floor(Xq); y0 floor(Yq); dx Xq - x0; dy Yq - y0; % 取四个邻域坐标越界时 clamp 到边界 x0c min(max(x0, 1), w); y0c min(max(y0, 1), h); x1c min(max(x0 1, 1), w); y1c min(max(y0 1, 1), h); idx00 sub2ind([h, w], y0c, x0c); idx01 sub2ind([h, w], y0c, x1c); idx10 sub2ind([h, w], y1c, x0c); idx11 sub2ind([h, w], y1c, x1c); I00 img(idx00); I01 img(idx01); I10 img(idx10); I11 img(idx11); out (1 - dx) .* (1 - dy) .* I00 ... dx .* (1 - dy) .* I01 ... (1 - dx) .* dy .* I10 ... dx .* dy .* I11; % 映射到原图范围之外的区域直接置黑 valid x0 1 x0 w - 1 y0 1 y0 h - 1; out(~valid) 0; end这里有个细节需要注意x0是向下取整得到的整数理论上x0的取值范围是[1, w-1]x01才是右下角的列号。但如果映射点刚好在图像右边界附近x0可能等于w那x01就超出了范围。我在代码里用了两次min(max())把越界坐标压缩到边界内保证sub2ind不会报错。这种处理方式相当于把边缘处的像素做了延伸实际效果比直接报错或者把整个边缘置黑更自然。双线性插值最大的变化是边缘变平滑了。斜线不再是生硬的台阶而是一段段有灰度过渡的斜边肉眼看上去舒服很多。代价是——边缘会被轻微模糊因为加权平均本质上是一个低通滤波操作它会损失一些高频细节。细纹理丰富的图像用双线性旋转后仔细对比能看出细节没有原来的锐利。2.4 双三次插值边缘与纹理的细节控双三次插值又叫三次卷积插值理论上它是三种算法里图像质量最高的。它利用了映射坐标周围4x4共16个邻域像素每个像素的权重由一个三次多项式核函数计算得到。常用的卷积核公式取a-0.5W(d) (a2)|d|³ - (a3)|d|² 1当 |d| ≤ 1 W(d) a|d|³ - 5a|d|² 8a|d| - 4a当 1 |d| 2 W(d) 0其他情况式中的d是映射点到某个邻域像素的带符号距离。|d|越小的像素权重越大距离超过2的像素完全不参与计算。这个核函数的曲线形状有点像加窗的sinc函数既能平滑过渡又能保留相对多的边缘信息所以主观画质比双线性好。function out bicubic_interp(img, Xq, Yq) % 双三次插值 [h, w] size(img); a -0.5; x0 floor(Xq); y0 floor(Yq); dx Xq - x0; dy Yq - y0; out zeros(size(Xq)); % 遍历 4x4 邻域偏移量 i 对应列方向j 对应行方向 for i -1:2 for j -1:2 xi x0 i; yi y0 j; % 越界 clamp 到边界 xi_c min(max(xi, 1), w); yi_c min(max(yi, 1), h); idx sub2ind([h, w], yi_c, xi_c); wx cubic_weight(dx - i, a); wy cubic_weight(dy - j, a); out out wx .* wy .* img(idx); end end % 越界区域置黑 valid x0 2 x0 w - 2 y0 2 y0 h - 2; out(~valid) 0; end function w cubic_weight(d, a) % 三次卷积核权重 d abs(d); w zeros(size(d)); m d 1; n d 1 d 2; w(m) (a 2) * d(m).^3 - (a 3) * d(m).^2 1; w(n) a * d(n).^3 - 5 * a * d(n).^2 8 * a * d(n) - 4 * a; end代码里的循环是两层4次循环方向别搞混i是列方向偏移j是行方向偏移。wx是用dx减i求出的这对应了“映射点与第i个邻域点在水平方向上的距离”wy同理对应垂直方向。双三次插值的运算量明显上来了每个输出像素要计算16个权重、16次乘加比双线性多了好几倍。但换来的是旋转后图像细节保留更好纹理边缘相对清晰没有双线性那种明显的模糊感。不过它也不是万能的——如果输入图本身质量不高或者旋转角度非常大导致大面积区域是插值生成而非原图真实信息双三次的感觉也不会神奇到哪去。3. 质量与性能对比实验3.1 测试图设计与实验方法只看原理和代码太抽象了。我用三张测试图做了对比实验一张是256×256的人工合成图含水平线、垂直线、斜线和同心圆一张是512×512的自然灰度图纹理细节比较丰富还有一张是带文字的截图用来观察边缘和文字的还原效果。合成图的生成代码很简单img zeros(256, 256); % 水平线和垂直线 img(128, :) 1; img(:, 128) 1; % 斜线45度 for k 1:256 img(k, 256 - k 1) 1; end % 同心圆 [X, Y] meshgrid(1:256, 1:256); r sqrt((X - 128.5).^2 (Y - 128.5).^2); img(abs(r - 40) 0.5) 1; img(abs(r - 80) 0.5) 1;这种图对旋转算法的观察特别友好直线够直斜线角度明确圆形弧度清晰任何锯齿、变形、模糊都逃不过眼睛。以30度旋转为例对三种算法逐个调用img_rgb imread(test.png); img_gray rgb2gray(img_rgb); r30_nearest my_rotate(img_gray, 30, nearest); r30_bilinear my_rotate(img_gray, 30, bilinear); r30_bicubic my_rotate(img_gray, 30, bicubic);然后通过imshow分别观察再截取图像中央的局部区域放大对比。3.2 三种算法的输出效果对照直接说肉眼能看到的差异。最近邻的结果最突出的就是锯齿。水平线和垂直线旋转30度后线条边缘呈现明显的阶梯状斜线更惨会出现忽粗忽细的断裂感同心圆在高曲率区域可以看到像素点拼成的多边形轮廓。整体观感偏“硬”线条像是用铅笔画在粗网格纸上。双线性插值的结果阶梯大幅减少线条边缘有了灰度过渡肉眼看起来平滑很多。但牺牲也明显线条整体变“糊”了。细看文字边缘会有淡出的灰边同心圆的弧线边缘不锐利。如果原图高频信息多这种模糊感会更突出。双三次插值的结果介于两者之间但更偏双线性一端——平滑度不错同时边缘的锐利度比双线性更好。文字轮廓更清晰直线边缘更干脆同心圆的弧形更接近理想圆。主观质量在三种算法里排第一。我做了个对比表方便直观参考评价维度最近邻双线性双三次边缘锯齿明显轻微最轻微细节保留区域值跳跃细节模糊细节较好线条平滑度差好最好文字可读性有断裂感略糊清晰计算复杂度最低中等最高需要注意双三次的“好”不是绝对的。对于锐利度已经被破坏的老照片、压缩率很高的JPEG图双三次和双线性的差距肉眼几乎看不出来但耗时差距却实打实摆在那。3.3 耗时与场景选型建议从算法复杂度推每个输出像素的计算量最近邻是一次比较和一次查表双线性要做8次浮点乘加双三次要做16个权重函数求值加16次乘加。实际跑起来双线性大约是最近邻的2到3倍耗时双三次大约是双线性的3到5倍比最近邻高一个数量级。我在普通配置的笔记本上用512×512灰度图旋转30度三种算法的耗时差距基本符合这个量级。双三次肉眼便能察觉的等待感最近邻则是瞬间出图。选型建议直接给结论工程预览、实时交互、嵌入式设备优先最近邻图快。通用处理、打印输出、日常图像操作选双线性性价比最高。高质量图像处理、图像超分前的预处理、视觉算法里需要最大限度保留纹理信息的场合选双三次。4. 实操中常见的坑与排查实录4.1 黑边、花斑和空洞我自己写第一版旋转代码时用的就是正向映射。结果旋转30度后图像里规律性地出现斜向黑色条纹像是被什么东西“啃”过一样。检查后发现这就是1.2节说的空洞问题——多个原图像素挤到同一个目标位置另一些目标位置没人管。解决方式很直接换成逆向映射扫一遍目标图像的每个像素它一定对应原图一个有效位置。改完之后空洞立刻消失。如果你在别人的代码里看到正向映射然后还有黑条纹的多半是没做“反向填充”或者只做了简单取整。判断一段旋转代码用的正向还是逆向最简单的方法是看遍历的是原图还是目标图遍历原图、往新图里填的是正向遍历新图、往原图里查的是逆向。4.2 “1像素偏移”之谜旋转中心如果直接用w/2、h/2出来的图整个会和正确结果差半个像素。这个偏移在视觉上不容易一眼看出来但如果你把旋转结果和imrotate的结果做差会看到图像边缘有半像素宽的灰度残影中心区域隐约有轮廓错位。原因前面提过MATLAB图像像素的存储位置是按“网格中心点”算的一行有w个像素第一个像素中心在1最后一个在w整行的对称中心在(w1)/2。所以无论原图还是目标图旋转中心都应该用(size1)/2。我在前几版代码里也吃过这个亏后来写了一个小函数验证把图像绕中心旋转360度再转回来如果中心坐标用的是(size1)/2旋转前后图像几乎完全重合用size/2就会有整体平移。这个方法推荐你也试一试。4.3 彩色图像扩展三通道前面的代码处理的是灰度图彩色图怎么转最简单的办法是把R、G、B三个通道拆开分别做一次旋转再合回去function rot_rgb my_rotate_rgb(img_rgb, angle_deg, method) R my_rotate(img_rgb(:, :, 1), angle_deg, method); G my_rotate(img_rgb(:, :, 2), angle_deg, method); B my_rotate(img_rgb(:, :, 3), angle_deg, method); rot_rgb cat(3, R, G, B); end反正三个通道共享同一套坐标映射各算各的像素值就行。如果想做点优化也可以先把RGB转成亮度-色度模型只对亮度通道做插值、色度通道用最近邻视觉差异很小但能省不少时间。这种技巧适合做实时彩色图像旋转时用不过前提是对颜色精度要求没那么苛刻。4.4 性能优化向量化是王道写MATLAB算法最忌讳的就是在像素级用双层for循环。我看到过不少版本代码用for row 1:rot_h、for col 1:rot_w逐像素计算坐标和插值一张1024×1024的图转30度要跑几十秒让人怀疑人生。正确做法就是把坐标和计算全部向量化。axis的meshgrid直接生成整个目标图的坐标矩阵坐标变换用矩阵点乘一次算完插值里sub2ind和权重公式都是矩阵运算全程不写一个像素级循环。除非是双三次那种4×4邻域循环那也只需要循环16次每次处理的还是整张图对应的矩阵不是像素级循环。还有一种优化是换精度。默认的double矩阵占8字节图像处理里很多环节其实用single就够了内存占用减半速度也会快一些。我在大图上旋转时会把图像转成single算完再转回uint8。4.5 与内置imrotate的结果差异写完手写版本你肯定会拿imrotate的结果做对比。发现不一样别慌。MATLAB自带的imrotate做了几个额外的处理默认填充方式是“边界外补0”但还支持crop保持原尺寸、loose自动扩边。内置函数用了更精细的坐标映射和边缘处理策略可能包含预滤波步骤。内置版本对边界像素的处理比手写版更讲究不会出现明显的边缘拉伸失真。所以手写算法的输出和imrotate有细微差异是正常的。只要你的视觉效果合理、没有严重伪影实验就算合格。真要做定量对比可以算旋转前后图像的峰值信噪比或者结构相似度指标差异在合理范围内即可。4.6 边界处理策略小结插值坐标落到图像边缘附近时怎么办我在三种插值函数里统一采用了“坐标clamp到边界”和“越界置黑”的组合策略。这里有几个可选的思路按效果排序补零填充越界部分当作0简单但会形成明显黑色背景视觉突兀。边缘复制越界坐标全取最近的边缘像素旋转后背景是原图边缘颜色看起来自然。镜像填充越界坐标按镜像反射到图像内部比边缘复制更平滑适合对边缘质量要求高的算法。课程实验用补零或者边缘复制足够了。要是做精细的图像拼接、纹理分析建议试试镜像填充边缘锯齿会小很多。我个人在写这个旋转函数时最深的体会是整个实现里最影响成败的不是插值公式而是坐标映射那一环。中心点算错、方向定义错、正逆映射搞反任何一个小问题都会让结果花掉。代码写完后多做对比实验、多画测试图比对着公式空想要有效得多。如果后面想继续扩展可以试试把旋转和缩放结合起来或者把变换矩阵推广成通用的仿射变换那又是一片新的玩法了。