matlab图像拼接的四种方法

 1、直接拼接，

 2、亮度调整后拼接，

 3、按距离比例融合，

 4、亮度调整后按距离比例融合



流程：

1。读入左，右图，并取出重合部分，并转化为亮度图

2。分别把每点的亮度值相加，得到一个比值

3。把比值 乘以 右图

4。再把左 各 右图 拼接

5。权重融合

左图------------重合区----------------------右图（相加）

100%-----------50%------------------------0%

0%--------------50%------------------------100%

 

%读入原图 （左 右）
img1=imread('2.jpg');
img2=imread('1.jpg');




 

 

%查找它们的SIFT特征,并返回匹配点对---------------------特征匹配 开始
[des1, des2] = siftMatch(img1, img2);
des1=[des1(:,2),des1(:,1)];%左右（x和y）交换 % 为过滤匹配准备参数
des2=[des2(:,2),des2(:,1)];%

%用 基础矩阵F 过滤匹配的特征点对
matchs = matchFSelect(des1, des2) %匹配位置索引（掩码）
des1=des1(matchs,:);%取出内点
des2=des2(matchs,:);

% 画出匹配特征点的连接线（好点）
drawLinedCorner(img1,des1,img2, des2) ;




 

1、直接拼接

 

[H,W,k]=size(img1);%图像大小
l_r=W-des1(1,2)+des2(1,2);%只取水平方向（第一个匹配点）重叠宽度


% 1、直接拼接-------------------------------------------------

%[H,W,k]=size(img1);
%l_r=405;%重叠宽度（W-宽 至 W）---如果不用特征匹配这里直接写重合区宽
L=W+1-l_r;%左边起点
R=W;%右边尾点
n=R-L+1;%重叠宽度：就是l_r
%直接拼接图
im=[img1,img2(:,n:W,:)];%1全图+2的后面部分
figure;imshow(im);title('直接拼接图');




 

2、亮度调整后拼接

 

%根据之前H矩阵找到的两幅图的重叠（l_r）部分
A=img1(:,L:R,:);
B=img2(:,1:n,:);
%A、B 是对应在两幅图中的重叠区域

% A=uint8(A); figure;imshow(A);
% B=uint8(B);figure;imshow(B);
% 
[ma,na,ka]=size(A);
I1=rgb2gray(A);%转换为灰度图像
I1=double(I1);%转换为双精度
v1=0;
I2= rgb2gray(B);
I2=double(I2);
v2=0;
for i=1:ma
    for j=1:na
        %I1(i,j)=0.59*A(i,j,1)+0.11*A(i,j,2)+0.3*A(i,j,3);%按点转化为灰度图
        v1=v1+I1(i,j);%所有亮度值相加（和）
        %I2(i,j)=0.59*B(i,j,1)+0.11*B(i,j,2)+0.3*B(i,j,3);
        v2=v2+I2(i,j);
    end
end

%figure;imshow(I1,[]);
%figure;imshow(I2,[]);

%亮度比例，并按比例调整第二个图
k=v1/v2;

BB2=img2(:,n:W,:)*k;%乘比值
im2=[img1,BB2];%拼接
figure;imshow(im2);title('调整亮度后拼接图');





 

3、按距离比例融合

 

%[H,Y,t]=size(im);
C=im;%继承前图
D=im2;%继承前图（亮度）
% n=拼缝宽;
%for i=1:H %少一重循环
    for j=1:n
        d=1-(j)/n;%disp(d);% 距离权重
        C(1:H,L+j,:)=d*A(1:H,j,:)+(1-d)*B(1:H,j,:);%互补融合
        D(1:H,L+j,:)=d*A(1:H,j,:)+(1-d)*B(1:H,j,:)*k;
    end
%end
C=uint8(C);
figure;imshow(C);title('直接融合拼接图');%3
D=uint8(D);
figure;imshow(D);title('亮度处理后融合拼接图');%4




 

4、亮度调整后按距离比例融合



完整的m文件：

 

%matlab图像拼接(四种方法)
% 1、直接拼接，
% 2、亮度调整后拼接，
% 3、按距离比例融合，
% 4、亮度调整后按距离比例融合

%流程：
%1。读入左，右图，并取出重合部分，并转化为亮度图
%2。分别把每点的亮度值相加，得到一个比值
%3。把比值 乘以 右图
%4。再把左 各 右图 拼接

clear;close all,clc;
%读入原图 （左 右）
img1=imread('2.jpg');
img2=imread('1.jpg');
% figure;imshow(img1);%显示
% figure;imshow(img2);

%查找它们的SIFT特征,并返回匹配点对---------------------特征匹配 开始
[des1, des2] = siftMatch(img1, img2);
des1=[des1(:,2),des1(:,1)];%左右（x和y）交换 为基础矩阵F 过滤匹配准备参数
des2=[des2(:,2),des2(:,1)];%

%用 基础矩阵F 过滤匹配的特征点对
matchs = matchFSelect(des1, des2) %匹配位置索引（掩码）
des1=des1(matchs,:);%取出内点
des2=des2(matchs,:);

% 画出匹配特征点的连接线（好点）
drawLinedCorner(img1,des1,img2, des2) ;
%------------------------------------------------------特征匹配 结束

[H,W,k]=size(img1);%图像大小
l_r=W-des1(1,2)+des2(1,2);%只取水平方向（第一个匹配点）重叠宽度


% 1、直接拼接-------------------------------------------------

%[H,W,k]=size(img1);
%l_r=405;%重叠宽度（W-宽 至 W）---如果不用特征匹配这里直接写重合区宽
L=W+1-l_r;%左边起点
R=W;%右边尾点
n=R-L+1;%重叠宽度：就是l_r
%直接拼接图
im=[img1,img2(:,n:W,:)];%1全图+2的后面部分
figure;imshow(im);title('直接拼接图');

% 2、亮度调整后拼接-------------------------------------------------
%根据之前H矩阵找到的两幅图的重叠（l_r）部分
A=img1(:,L:R,:);
B=img2(:,1:n,:);
%A、B 是对应在两幅图中的重叠区域

% A=uint8(A); figure;imshow(A);
% B=uint8(B);figure;imshow(B);
% 
[ma,na,ka]=size(A);
I1=rgb2gray(A);%转换为灰度图像
I1=double(I1);%转换为双精度
v1=0;
I2= rgb2gray(B);
I2=double(I2);
v2=0;
for i=1:ma
    for j=1:na
        %I1(i,j)=0.59*A(i,j,1)+0.11*A(i,j,2)+0.3*A(i,j,3);%按点转化为灰度图
        v1=v1+I1(i,j);%所有亮度值相加（和）
        %I2(i,j)=0.59*B(i,j,1)+0.11*B(i,j,2)+0.3*B(i,j,3);
        v2=v2+I2(i,j);
    end
end

%figure;imshow(I1,[]);
%figure;imshow(I2,[]);

%亮度比例，并按比例调整第二个图
k=v1/v2;

BB2=img2(:,n:W,:)*k;%乘比值
im2=[img1,BB2];%拼接
figure;imshow(im2);title('调整亮度后拼接图');


% 3、按距离比例融合-------------------------------------------------
% 4、亮度调整后按距离比例融合----------------------------------------

% 图像融合消除拼接缝隙

%用的渐入渐出融合即：距离权重融合

%[H,Y,t]=size(im);
C=im;%继承前图
D=im2;%继承前图（亮度）
% n=拼缝宽;
%for i=1:H %少一重循环
    for j=1:n
        d=1-(j)/n;%disp(d);% 距离权重
        C(1:H,L+j,:)=d*A(1:H,j,:)+(1-d)*B(1:H,j,:);%互补融合
        D(1:H,L+j,:)=d*A(1:H,j,:)+(1-d)*B(1:H,j,:)*k;
    end
%end
C=uint8(C);
figure;imshow(C);title('直接融合拼接图');%3
D=uint8(D);
figure;imshow(D);title('亮度处理后融合拼接图');%4


 



特征点匹配函数取自一个 siftDemoV4 包

 

 

matlab 7.0 运行

感谢 浮华三生 的评论！

我把siftMatch放在后面，并去掉matchFSelect函数，用上包上的findHomography函数。

先在 这里下 http://www.cs.ubc.ca/~lowe/keypoints/siftDemoV4.zip

%查找它们的SIFT特征,并返回匹配点对---------------------特征匹配 开始
[des1, des2] = siftMatch(img1, img2);

%为单应矩阵过滤匹配准备参数
pts1=des1';pts2=des2';

%单应矩阵过滤匹配
[Ht matchs] = findHomography(pts1,pts2);%如无该函数，自行搜索“imMosaic”或“findHomography imMosaic”下载
pts1=pts1(:,matchs);%取出内点
pts2=pts2(:,matchs);
des1=pts1';%格式转回
des2=pts2';

% 画出匹配特征点的连接线（好点）
drawLinedCorner(img1,des1,img2, des2) ;
%------------------------------------------------------特征匹配 结束

替换前面部分

修改match函数为siftMatch:

% 此函数读取两幅图像，查找它们的SIFT特征
% 仅当匹配的距离小于与第二个最接近匹配的距离的阈值时，才接受匹配。
% 它返回两个图像的匹配点，matchLoc1 = [x1,y1;x2,y2;...]
%
% 

function [matchLoc1 matchLoc2] = siftMatch(img1, img2)

% 在每个图像查找 SIFT 特征点
[des1, loc1] = sift(img1);
[des2, loc2] = sift(img2);

% 对于MATLAB中的效率,计算单位向量之间的点乘比欧式距离更快捷。请注意: 
% 注意角的比率(单位矢量点积的反余弦)是小角度的欧氏距离之比的近似值。
%
% distRatio: 在这两队匹配中只保留矢量角从最近的第二近邻的比值小于distRatio的。
distRatio = 0.7;   

% 在第一个图像中的每个描述符，选择它的匹配到第二个图像。
des2t = des2';                          % 预计算矩阵转置 
matchTable = zeros(1,size(des1,1));
for i = 1 : size(des1,1)
   dotprods = des1(i,:) * des2t;        % 计算点积向量
   [vals,indx] = sort(acos(dotprods));  % 取逆余弦和排序结果

   % 看看和最近的邻居角比率小于distRatio。
   if (vals(1) < distRatio * vals(2))
      matchTable(i) = indx(1);
   else
      matchTable(i) = 0;
   end
end
% 保存匹配数据表

num = sum(matchTable > 0);
fprintf('找到 %d 对匹配点.\n', num);

idx1 = find(matchTable);
idx2 = matchTable(idx1);
x1 = loc1(idx1,2);
x2 = loc2(idx2,2);
y1 = loc1(idx1,1);
y2 = loc2(idx2,1);

matchLoc1 = [y1,x1];%把y坐标放前面
matchLoc2 = [y2,x2];

end

说明一下：

http://www.cs.ubc.ca/~lowe/keypoints/siftDemoV4.zip 还有findHomography的 https://github.com/erogol/Image_Stiching

 

 

 

 

说到图像处理，第一步就是图像读取。Matlab最简单的就是imread函数，本节介绍imread的用法以及容易出错的地方


正如上图所示，在Matlab文档中，imread包括以上几种用法，但并不要求都要掌握，个人觉得会用一两种，其他了解语法即可。

我们来介绍下最最常用的语句A = imread(filename)

我们来读取一张图片
>> a = imread('凝视.jpg');
>> imtool(a)



如图所示，首先注意的是语法正确A = imread(‘凝视.jpg’);
一.正确示范。

他的意思是在当前路径里把文件名为“凝视.jpg”的图片数据读取到A中保存，那我们可以看到在图片最右边Workspace区域放着一些数据，这就是A的数据，我们看到这张图片是3405933大小，意思是340行，593列，3通道（RGB）的图片，右边的UINT8表示8位无符号的整型类型。（后面的imtool语句用于显示图片，这个将在后续细谈）

补充一点小知识：

想要清空Command Window,输入命令clc

想要清空Workspace,输入命令clear

想要关闭所有打开的窗口，输入命令close all

想要观看图像信息，用whos
二.踩坑

好了，既然我们知道正确的写法，那我们来试试有什么坑。（敢于试错是干这行的一个优秀品质）

1.为什么要分号？

因为matlab是按行编译，一行一行编译，不写分好，就会直接出来编译结果，给你们看个例子就明白。


上面创建一个a和b矩阵，a矩阵不用分号结尾，窗口直接显示内容，而b矩阵用了分号，并不显示内容，但是可以看到编译后，在workspace已经创建两个数组矩阵，我们也可以看到，点击变量名可以看到最上面有具体的数据。同理，如果我们读取图片a = imread(‘凝望.jpg’)不写分号，那么窗口就会出现一大波数据，刷刷刷地跳出来，那酸爽，那些数据就是保存在数组里面的像素了。

小知识：

如果你不写变量名，比如>>imread(‘凝望.jpg’);它很默认地给你个名字叫：ans

当你想重新写一句语句和上面相同或相似的话，可以按键盘的上箭头按钮


这个快捷方式可以帮助你快速地修改语句，好用。
2.英文半角符号

这个学过编程的都应该知道吧，别整个中文“。”，也千万别用英文的全角符号“．”,至于为什么，我不知道。Matlab会出现：File “凝视．jpg” does not exist.
3.路径

这个是常会犯的，你记得你有那图片，也记得图片名称，但是你没有把他放在当前的路径里你叫人家怎么找，Matlab还没有那么强大到可以全搜索你电脑里面的图片文件。同样编译会出现 ： does not exist。但是还有补救的方法，你可以给它指明可以去哪里找，例如:


我把图片‘凝视.jpg’放在D盘的build文件下面，结果不断地出错出错如图所示，乍一看，才发现把凝视写成凝望了，再一看，把build写成bulid，哈哈，我也只是个初学者，稍粗心就犯错，借此也说明写代码确实得聚精会神啊，这点BUG还容易改，但是如果你做的是大项目，写成百上千行代码，因为拼写错误而花了几个小时那实在是亏大了。
好了，本次文件读写就写到这里，有什么疑问可以评论大家一起讨论学习，也许我们会有思想的火花。可能说的都是很琐碎很简单的东西，只要能给你一点点收获，这篇博客就有价值了。下节继续介绍其他函数。感谢观看

基于MATLAB图像处理的中值滤波、均值滤波以及高斯滤波的实现与对比

1.  背景知识

中值滤波法是一种非线性平滑技术，它将每一像素点的灰度值设置为该点某邻域窗口内的所有像素点灰度值的中值.

中值滤波是基于排序统计理论的一种能有效抑制噪声的非线性信号处理技术，中值滤波的基本原理是把数字图像或数字序列中一点的值用该点的一个邻域中各点值的中值代替，让周围的像素值接近的真实值，从而消除孤立的噪声点。

方法是用某种结构的二维滑动模板，将板内像素按照像素值的大小进行排序，生成单调上升（或下降）的为二维数据序列。二维中值滤波输出为g（x,y）=med{f(x-k,y-l),(k,l∈W)} ，其中，f(x,y)，g(x,y)分别为原始图像和处理后图像。W为二维模板，通常为3*3，5*5区域，也可以是不同的的形状，如线状，圆形，十字形，圆环形等。

 
			
 
			
 
			
 
			
 
		
 
			

			
2
			
			

			
4
			
			

			
8
			
			
 
		
 
			

			
1
			
			

			
3
			
			

			
9
			
			
 
		
 
			

			
5
			
			

			
7
			
			

			
6
			
			
 
		
 
			
 
			
 
			
 
			
 
		
g（x,y）=med{f(x-k,y-l),(k,l∈W)}

g =med[2,4,8;1,3,9;5,7,6] = 5

中值滤波后的结果：

 
			
 
			
 
			
 
			
 
		
 
			
 
			
 
			
 
			
 
		
 
			
 
			

			
5
			
			
 
			
 
		
 
			
 
			
 
			
 
			
 
		
 
			
 
			
 
			
 
			
 
		
均值滤波是典型的线性滤波算法，它是指在图像上对目标像素给一个模板，该模板包括了其周围的临近像素（以目标像素为中心的周围8个像素，构成一个滤波模板，即去掉目标像素本身），再用模板中的全体像素的平均值来代替原来像素值。

均值滤波也称为线性滤波，其采用的主要方法为邻域平均法。线性滤波的基本原理是用均值代替原图像中的各个像素值，即对待处理的当前像素点（x，y），选择一个模板，该模板由其近邻的若干像素组成，求模板中所有像素的均值，再把该均值赋予当前像素点（x，y），作为处理后图像在该点上的灰度g（x，y），即g（x，y）=1/m ∑f（x，y） m为该模板中包含当前像素在内的像素总个数。

均值滤波本身存在着固有的缺陷，即它不能很好地保护图像细节，在图像去噪的同时也破坏了图像的细节部分，从而使图像变得模糊，不能很好地去除噪声点。

 
			
 
			
 
			
 
			
 
		
 
			

			
2
			
			

			
4
			
			

			
8
			
			
 
		
 
			

			
1
			
			

			
3
			
			

			
9
			
			
 
		
 
			

			
5
			
			

			
7
			
			

			
6
			
			
 
		
 
			
 
			
 
			
 
			
 
		
g（x，y）=1/m ∑f（x，y）

g = (1/8)*(2+4+8+1+9+5+7+6)= 5

均值滤波后的结果：

 
			
 
			
 
			
 
			
 
		
 
			
 
			
 
			
 
			
 
		
 
			
 
			

			
5
			
			
 
			
 
		
 
			
 
			
 
			
 
			
 
		
 
			
 
			
 
			
 
			
 
		
 

高斯滤波是一种线性平滑滤波，适用于消除高斯噪声，广泛应用于图像处理的减噪过程。通俗的讲，高斯滤波就是对整幅图像进行加权平均的过程，每一个像素点的值，都由其本身和邻域内的其他像素值经过加权平均后得到。高斯滤波的具体操作是：用一个模板（或称卷积、掩模）扫描图像中的每一个像素，用模板确定的邻域内像素的加权平均灰度值去替代模板中心像素点的值。

2.  MATLAB实现

源码：

%%-------------------------------------------------------------------------------------------

%% 2018/01/03

%% lee

%% 137194782@qq.com

%% 微信公众号：FPGA开源工作室

%%--------------------------------------------------------------------------------------------

clear all;

clc;

M = imread('timg.jpg');         %读取MATLAB中的名为timg的图像   

figure,imshow(M);                %显示原始图像

title('original');

gray = rgb2gray(M);

figure,imshow(gray);                         %显示灰度图像

title('gray');

P1 = imnoise(gray,'gaussian',0.02);     %加入高斯躁声 

figure,imshow(P1);                        %加入高斯躁声后显示图像

title('gaussiannoise');

P2 = imnoise(gray,'salt& pepper',0.02); %加入椒盐躁声

figure,imshow(P2);                        %加入椒盐躁声后显示图像  

title('salt& pepper noise');

g = medfilt2(P1);                       %对高斯躁声中值滤波

figure,imshow(g);

title('medfiltergaussian');

h = medfilt2(P2);                       %对椒盐躁声中值滤波

figure,imshow(h);

title('medfiltersalt & pepper noise');

a=[1 1 1                               %对高斯躁声算术均值滤波

   1 11

   1 11];

l=1/9*a;                          

k = conv2(double(P1),double(l));                        

figure,imshow(k,[]);

title('arithmeticfiltergaussian');                                          

d = conv2(double(P2),double(l));           %对椒盐躁声算术均值滤波              

figure,imshow(d,[]);

title('arithmeticfiltersalt & pepper noise');

sigma=8;%标准差大小  

window=double(uint8(3*sigma)*2+1);%窗口大小一半为3*sigma  

H=fspecial('gaussian',window, sigma);%fspecial('gaussian',hsize, sigma)产生滤波模板   

img_gauss=imfilter(P1,H,'replicate'); %为了不出现黑边，使用参数'replicate'（输入图像的外部边界通过复制内部边界的值来扩展）  

figure, imshow(img_gauss);

title('gaussianfilting gauss noise');

img_salt=imfilter(P2,H,'replicate'); 

figure, imshow(img_salt);

title('gaussianfilting salt pepper noise');

 

结果展示：

 



 

 

原始图像

 

 



 

灰度图像

 



 

 

加入高斯噪声的灰度图像



 

 

 

加入椒盐噪声的灰度图像

 



 

 

经过中值滤波后的高斯噪声灰度图像

 



 

经过中值滤波后的椒盐噪声灰度图像

 



 

经过均值滤波的高斯噪声灰度图像

 



 

经过均值滤波的椒盐噪声灰度图像

 



 

经过高斯滤波的高斯噪声灰度图像



经过高斯滤波的椒盐噪声的灰度图像

结果分析：图像经过中值滤波后，高斯噪声没有被完全去除，椒盐噪声几乎被完全去除效果较好。经过均值滤波后不管是高斯噪声还是椒盐噪声大部分都没有被去除，只是稍微模糊化。经过高斯滤波后，高斯噪声和椒盐噪声几乎被很大程度的模糊化，原图好像被加上了一层蒙版。

 

欢迎大家关注我的微信公众号FPGA开源工作室

FPGA自习学院将不断更新和总结FPGA相关学习资料，书籍，实例工程和视频。

欢迎大家加入FPGA自习学院，一起学习一起成长。



 

 

 

 

 

                                                           Matlab图像旋转

####以坐标原点为中心旋转的原理：

点p0p_0p0​绕坐标原点逆时针方向旋转$\theta$角度得到点p1p_1p1​.


从①②可以得到旋转变换的矩阵公式：
[x1y11]\begin{bmatrix}x_1&y_1&1\end{bmatrix}[x1​​y1​​1​]=[x0y01]\begin{bmatrix}x_0&y_0&1\end{bmatrix}[x0​​y0​​1​][cos(θ)sin(θ)0−sin(θ)cos(θ)0001]\begin{bmatrix}cos(\theta)&sin(\theta)&0\\-sin(\theta)&cos(\theta)&0\\0&0&1\end{bmatrix}⎣⎡​cos(θ)−sin(θ)0​sin(θ)cos(θ)0​001​⎦⎤​
还可以得到逆运算矩阵公式：
[x0y01]\begin{bmatrix}x_0&y_0&1\end{bmatrix}[x0​​y0​​1​]=[x1y11]\begin{bmatrix}x_1&y_1&1\end{bmatrix}[x1​​y1​​1​][cos(θ)−sin(θ)0sin(θ)cos(θ)0001]\begin{bmatrix}cos(\theta)&-sin(\theta)&0\\sin(\theta)&cos(\theta)&0\\0&0&1\end{bmatrix}⎣⎡​cos(θ)sin(θ)0​−sin(θ)cos(θ)0​001​⎦⎤​

####以任意图形中心点为坐标原点旋转原理：

从上图可知以任意图形中心点为坐标原点旋转我们需要三步：

(1)将坐标系Ⅰ变成坐标系Ⅱ

(2)在坐标系Ⅱ中旋转$\theta$角

(3)将坐标系Ⅱ变成坐标系Ⅰ
#####(1)将坐标系Ⅰ变成坐标系Ⅱ

由Figure1得到Figure2可知，变换矩阵为：
[xⅡyⅡ1]\begin{bmatrix}xⅡ\\yⅡ\\1\end{bmatrix}⎣⎡​xⅡyⅡ1​⎦⎤​=[xⅠyⅠ1]\begin{bmatrix}xⅠ\\yⅠ\\1\end{bmatrix}⎣⎡​xⅠyⅠ1​⎦⎤​[1000−10−0.5w0.5h1]\begin{bmatrix}1&0&0\\0&-1&0\\-0.5w&0.5h&1\end{bmatrix}⎣⎡​10−0.5w​0−10.5h​001​⎦⎤​=[xⅠ−0.5w−(yⅠ−0.5h)1]\begin{bmatrix}xⅠ-0.5w\\-(yⅠ-0.5h)\\1\end{bmatrix}⎣⎡​xⅠ−0.5w−(yⅠ−0.5h)1​⎦⎤​
#####(2)在坐标系Ⅱ中旋转$\theta$角

见上面以坐标原点为中心旋转的原理
#####(3)将坐标系Ⅱ变成坐标系Ⅰ

由Figure3得到Figure4可知，变换矩阵为(其实就是(1)中变换矩阵的逆变换)：
[xⅠyⅠ1]\begin{bmatrix}xⅠ\\yⅠ\\1\end{bmatrix}⎣⎡​xⅠyⅠ1​⎦⎤​=[xⅡyⅡ1]\begin{bmatrix}xⅡ\\yⅡ\\1\end{bmatrix}⎣⎡​xⅡyⅡ1​⎦⎤​[1000−100.5nW0.5nH1]\begin{bmatrix}1&0&0\\0&-1&0\\0.5nW&0.5nH&1\end{bmatrix}⎣⎡​100.5nW​0−10.5nH​001​⎦⎤​=[xⅡ+0.5nW−yⅡ+0.5nH1]\begin{bmatrix}xⅡ+0.5nW\\-yⅡ+0.5nH\\1\end{bmatrix}⎣⎡​xⅡ+0.5nW−yⅡ+0.5nH1​⎦⎤​
从而将三步合起来，可以得到以任意图像中心为坐标原点旋转的变换矩阵：
[x1y11]\begin{bmatrix}x_1\\y_1\\1\end{bmatrix}⎣⎡​x1​y1​1​⎦⎤​=[x0y01]\begin{bmatrix}x_0\\y_0\\1\end{bmatrix}⎣⎡​x0​y0​1​⎦⎤​[1000−10−0.5w0.5h1]\begin{bmatrix}1&0&0\\0&-1&0\\-0.5w&0.5h&1\end{bmatrix}⎣⎡​10−0.5w​0−10.5h​001​⎦⎤​[cos(θ)sin(θ)0−sin(θ)cos(θ)0001]\begin{bmatrix}cos(\theta)&sin(\theta)&0\\-sin(\theta)&cos(\theta)&0\\0&0&1\end{bmatrix}⎣⎡​cos(θ)−sin(θ)0​sin(θ)cos(θ)0​001​⎦⎤​
[1000−100.5nW0.5nH1]\begin{bmatrix}1&0&0\\0&-1&0\\0.5nW&0.5nH&1\end{bmatrix}⎣⎡​100.5nW​0−10.5nH​001​⎦⎤​=[x0y01]\begin{bmatrix}x_0\\y_0\\1\end{bmatrix}⎣⎡​x0​y0​1​⎦⎤​[cos(θ)−sin(θ)0sin(θ)cos(θ)0−0.5wcos(θ)−0.5hsin(θ)+0.5nW−0.5wsin(θ)+0.5hcos(θ)+0.5nH1]\begin{bmatrix}cos(\theta)&-sin(\theta)&0\\sin(\theta)&cos(\theta)&0\\-0.5wcos(\theta)-0.5hsin(\theta)+0.5nW&-0.5wsin(\theta)+0.5hcos(\theta)+0.5nH&1\end{bmatrix}⎣⎡​cos(θ)sin(θ)−0.5wcos(θ)−0.5hsin(θ)+0.5nW​−sin(θ)cos(θ)−0.5wsin(θ)+0.5hcos(θ)+0.5nH​001​⎦⎤​
它的逆变换为：
[x0y01]\begin{bmatrix}x_0\\y_0\\1\end{bmatrix}⎣⎡​x0​y0​1​⎦⎤​=[x1y11]\begin{bmatrix}x_1\\y_1\\1\end{bmatrix}⎣⎡​x1​y1​1​⎦⎤​[1000−10−0.5nW0.5nH1]\begin{bmatrix}1&0&0\\0&-1&0\\-0.5nW&0.5nH&1\end{bmatrix}⎣⎡​10−0.5nW​0−10.5nH​001​⎦⎤​[cos(θ)−sin(θ)0sin(θ)cos(θ)0001]\begin{bmatrix}cos(\theta)&-sin(\theta)&0\\sin(\theta)&cos(\theta)&0\\0&0&1\end{bmatrix}⎣⎡​cos(θ)sin(θ)0​−sin(θ)cos(θ)0​001​⎦⎤​
[1000−100.5w0.5h1]\begin{bmatrix}1&0&0\\0&-1&0\\0.5w&0.5h&1\end{bmatrix}⎣⎡​100.5w​0−10.5h​001​⎦⎤​=[x0y01]\begin{bmatrix}x_0\\y_0\\1\end{bmatrix}⎣⎡​x0​y0​1​⎦⎤​[cos(θ)sin(θ)0−sin(θ)cos(θ)0−0.5nWcos(θ)+0.5nHsin(θ)+0.5w−0.5nWsin(θ)−0.5nHcos(θ)+0.5h1]\begin{bmatrix}cos(\theta)&sin(\theta)&0\\-sin(\theta)&cos(\theta)&0\\-0.5nWcos(\theta)+0.5nHsin(\theta)+0.5w&-0.5nWsin(\theta)-0.5nHcos(\theta)+0.5h&1\end{bmatrix}⎣⎡​cos(θ)−sin(θ)−0.5nWcos(θ)+0.5nHsin(θ)+0.5w​sin(θ)cos(θ)−0.5nWsin(θ)−0.5nHcos(θ)+0.5h​001​⎦⎤​
####Matlab编程实现

#####(1)Matlab自带函数实现图像任意角度旋转

######旋转函数介绍：

B=imrotate(A,angle,method, ‘crop’)

  angle   :旋转角度，单位为度，角度为正值时逆时针旋转

method  :该可选参数为imrotate()指定插值方法

‘crop’:裁剪旋转后增大的图像
clc                                 
I=imread('potted-plantsk.jpg');
figure,imshow(I);
title('srcImage');
I1=imrotate(I,30);                  %旋转30°
I2=imrotate(I,30,'crop');           %旋转30°，并剪切图像，使得到的图像和原图像大小一致
I3=imrotate(I,30,'bilinear','crop');%双线性插值法旋转30°，并剪切图像，使得到的图像和原图像大小一致
figure,imshow(I1);
title('I1');
figure,imshow(I2);
title('I2');
figure,imshow(I3);
title('I3');

程序运行结果如下：

#####(2)自编myimrotate()函数实现图像任意角度旋转
function [ A ] = myimrotate(B,degree)                                 %定义旋转函数，degree为要旋转的角度
[r,c,d]=size(B);                                                      %获取输入图像B的行r、列c和通道数d,为了旋转彩色图像所以有必要得到通道数d
nH=round(r*abs(cosd(degree))+c*abs(sind(degree)));                    %旋转图像后得到的新高度，“round()函数四舍五入“
nW=round(c*abs(cosd(degree))+r*abs(sind(degree)));                    %旋转图像后得到的新宽度
A=zeros(nH,nW,d);                                                     %定义生成目标图像的行列以及通道数
M1=[1 0 0;0 -1 0;-0.5*nW 0.5*nH 1 ];                                  %坐标系变换矩阵M1
M2=[cosd(degree) -sind(degree) 0;sind(degree) cosd(degree) 0;0 0 1];  %角度旋转变换矩阵M2，我用的是顺时针方向
M3=[1 0 0;0 -1 0;0.5*c 0.5*r 1];                                      %坐标系变换矩阵M3
    for i=1:nW
        for j=1:nH
            temp=[i j 1]*M1*M2*M3;                                    %得到旋转后的矩阵temp
            y=temp(1,2);                                              %y取矩阵temp的第一行第二列,y对应j，为高度
            x=temp(1,1);                                              %x取矩阵temp的第一行第一列,x对应i，为宽度
            y=round(y);                                               %y四舍五入取整
            x=round(x);                                               %x四舍五入取整
           if(x>=1&&x<=c)&&(y>=1&&y<=r)                               %判断的得到的(x,y)点是否在原图像上
               A(j,i,:)=B(y,x,:);                                     %将原图像的像素点赋值给对应的旋转后图像上的点
           end                                                        %（”有人疑惑为啥不是A(i,j,:)=B(x,y,:);因为i,x对应的是列，即宽，而j,y对应的是行，即高“），我这里以x为横坐标，y为竖向纵坐标
        end
    end
end

调用函数：
clc                                 
I=imread('potted-plantsk.jpg');
figure,imshow(I);
title('srcImage');
I1=myimrotate(I,30);     %调用myimrotate()函数旋转30° 
I2=myimrotate(I,-90);     %调用myimrotate()函数旋转-90°
figure,imshow(uint8(I1));
title('旋转30°：I1');
figure,imshow(uint8(I2));
title('旋转-90°：I2');

程序运行结果：

补充一下：

sin/cos/tan                     三角函数，角度单位为弧度

asin/acos/atan                反三角函数，角度单位为弧度
sind/cosd/tand               三角函数，角度单位为度

asind/acosd/atand          反三角函数，角度单位为度
补充：坐标系平移推导