MATLAB图像拼接技术研究

在实际应用中，经常会用到超过人眼视野范围甚至是全方位的高分辨率图像，普通数码相机的视野范围往往难以满足要求。为了得到大视野范围的图像，人们使用广角镜头和扫描式相机进行拍摄。但这些设备往往价格昂贵、使用复杂，而且在一幅低分辨率的图像中得到超宽视角会损失景物中物体的分辨率，所以，使用广角镜头和扫描式相机只能部分的解决这一问题。此外，广角镜头的图像边缘会难以避免的产生扭曲变形，不利于一些场合的应用。为了在不降低图像分辨率的条件下获取大视野范围的图像，人们提出了图像拼接技术，将普通图像或视频图像进行无缝拼接，得到超宽视角甚至360度的全景图，这样就可以用普通数码相机实现场面宏大的景物拍摄。利用计算机进行匹配，将多幅具有重叠关系的图像拼合成为一幅具有更大视野范围的图像，这就是图像拼接的目的。

本文的研究内容与组织结构

本文研究的重点是使用数码相机或手机手持拍摄的图像的全自动拼接。算法的基本流程如下：

1. 读取n副连续有重叠部分的图像，在n副图像中检测SIFT特征，并用SIFT

特征描述子对其进行描述。

1. 匹配相邻图像的特征点，并根据特征点向量消除误匹配。
2. 使用RANSAC方法，确定变换参数。
3. 图像融合

简介

如何确定矩阵M？ 首先确定一定数量的特征点，利用特征点的匹配给出图像变换的估计初值，最后通过递归算法得到最后的变换。

特征提取与匹配

采用SIFT算法进行特征点的提取与匹配。主要包含4个步骤：

1. 建立尺度空间,寻找候选点
2. 精确确定关键点,剔除不稳定点
3. 确定关键点的方向
4. 提取特征描述符

本文中matlab实现代码如下:

im11 = img{ccc-1};

im22 = img{ccc};

[im1, des1, loc1] = sift2(im11);

[im2, des2, loc2] = sift2(im22);

distRatio = ;

des2t = des2′;

match(size(des1,1)) = 0;

for i = 1 : size(des1,1)

dotprods = des1(i,:) * des2t;

[vals,indx] = sort(acos(dotprods));

if (vals(1) < distRatio * vals(2))

match(i) = indx(1);

end

end

im3 = appendimages(im1,im2);

figure(‘Position’, [100 100 size(im3,2) size(im3,1)]);

colormap(‘gray’);

imagesc(im3);

hold on;

cols1 = size(im1,2);

hold off;

num = sum(match > 0);

fprintf(‘Found %d matches.\n’, num);

其中sift2封装了sift特征点提取，返回值des是特征描述子；appendimages函数实现了两幅图像边挨边拼接。详细代码见和

消除误匹配

假设同一个相机拍出的图像尺度差别比较小。因此对于两对匹配点（P1，Q1）

（P2，Q2），向量P1P2, Q1Q2的夹角比较小，同时|P1P2| 和 |Q1Q2| 的值也比较接近。我们记录一个匹配点与其他匹配点形成的向量，并记录待匹配图像中的对应匹配点与其他相应匹配点形成的向量，通过比较，记录该匹配点通过匹配阈值的个数，如果超过匹配点总数的4/5，则该匹配点是正确匹配。否则为误匹配。

图1，消除误匹配

matlab实现代码如下:

KeyPoints(num,4) = 0;

j = 1;

for i = 1: size(des1,1)

if (match(i) > 0)

KeyPoints(j,1) = loc1(i,2);

KeyPoints(j,2) = loc1(i,1);

KeyPoints(j,3) = loc2(match(i),2);

KeyPoints(j,4) = loc2(match(i),1);

j = j+1;

end

end

yuzhi_y = min(size(im11,1),size(im22,1));

yuzhi_x = min(size(im11,2),size(im22,2));

yuzhi = sqrt(yuzhi_y^2+ yuzhi_x^2);

PointDis(num) = 0;

for i = 1:size(KeyPoints,1)

for j = 1:size(KeyPoints,1)

a = [abs(KeyPoints(i,1) – KeyPoints(j,1)),abs(KeyPoints(i,2) – KeyPoints(j,2))];

b = [abs(KeyPoints(i,3) – KeyPoints(j,3)),abs(KeyPoints(i,4) – KeyPoints(j,4))];

xiangsi = dot(a,b)/(norm(a)*norm(b));

if (xiangsi<cos(pi/4)||norm(a) – norm(b) >yuzhi/10)

PointDis(i) = PointDis(i) +1;

end

end

end

j = 0;

for i=1:num

if (PointDis(i) >num/6 )

KeyPoints(i,1) = inf;

j=j+1;

end

end

KeyPoints2(num -j,4) = 0;

z = 1;

for i=1:num

if (KeyPoints(i,1)~=inf)

KeyPoints2(z,:) = KeyPoints(i,:);

z= z+1;

end

end

KeyPoints = KeyPoints2;

num = size(KeyPoints,1);

for i = 1: num

line([KeyPoints(i,1) KeyPoints(i,3)+cols1], …

[KeyPoints(i,2) KeyPoints(i,4)], ‘Color’, ‘c’);

end

确定变换参数

我们用平面透视变换来表示图像间的几何关系：

为了恢复M中的8个参数，需要4对匹配点来求解。我们使用RANSAC(Random Sample Consensus)方法来寻找两幅图像间最佳的4对匹配点，得到相应的H阵，并剔除误配点。

前面自动提取和匹配得到的初始匹配特征点对集合中难免存在误配点，此外特征点在提取时也存在一定的精度误差，因此需要一种容错能力很强的算法来过滤初始匹配特征点对集合，而RANSAC正是这样的一种算法，其具体内容可描述如下：

给定N个数据点组成的数据集合P，假设这N个数据点中的绝大部分是由一个参数未知的特定模型产生，而该模型至少需要n个数据点来求解，且N>n，则可以通过下面的迭代过程求解该模型的参数。

将以下步骤运行k次： ’

(1)从P中随机选取n个数据点的子集S。

(2)由选取的这n个数据点计算出一个模型M。

(3)对数据集合中其余的N一n个数据点，计算它们与模型肘之间的距离，记录P中在M的某个误差允许范围内的数据点的个数c。

迭代(1)～(3)k次后，对应最大c值的模型即为所求，数据集合P中的这c个数据即为内点，其余的N—C个数据即为外点。

matlab实现代码如下:

DDmin =5000;

for ll = 1:6

zNum = 1;

while (1)

if (zNum>num*6)

break;

end

zNum = zNum+1;

DD(15) = 0;

H(3,3,15) = 0;

for k = 1:15

while (1)

suiji = round(rand(1,4)*(num-1) +1);

if (sum(find(suiji == inf)) == 0)

break;

end

end

j =1;

for i = 1:4

x2 = KeyPoints(suiji(i),1); y2 = KeyPoints(suiji(i),2);

A(j,:) = [x , y,1,0,0,0,-x*x2, -x2*y];

D(j) = x2;

j= j +1;

A(j,:) = [0,0,0,x,y,1,-x*y2,-y2*y];

D(j) = y2;

j = j+1;

%B(i-1,:)=[x,y,x2,y2];

end

y = inv(A)*D’;

if (sum(H == inf)>0)

continue;

else

for in1 = 1:num

TT = H(:,:,k)*[KeyPoints(in1,3); KeyPoints(in1,4); 1];

DD(k) = DD(k)+sqrt((TT(1) – KeyPoints(in1,1))^2+(TT(2) – KeyPoints(in1,2))^2);

end

end

end

[val3,indx3] = sort(DD);

i = 1;

while (1)

if (H(1,1,indx3(i))==inf||DD(indx3(i))==0)

i = i + 1;

else

break;

end

end

if (DD(indx3(i)) <=DDmin)

DDmin =DD(indx3(i));

Hfinal = H(:,:,indx3(i));

break;

end

end

end

for tc = 1:ccc-1

if (ccc-tc-1==0)

break;

else

Hfinal = HfinalMax(:,:,ccc-tc-1)*Hfinal;

end

end

图像融合

1.图像旋转之后可能新的像素点并不在实际的像素位置，因此需要插值处理。建议使用双线性插值。

2.对于重合的部分的像素值可以有如下计算办法：

图2 图像融合方法

我们采用method1

matlab实现代码如下:

edgept1 = round(Hfinal*[size(im22,2);size(im22,1);1]);

edgept2 = round(Hfinal*[1;size(im22,1);1]);

edgept3 = round(Hfinal*[1;1;1]);

edgept4 = round(Hfinal*[size(im22,2);1;1]);

height = max([size(im11,1),edgept1(2),edgept2(2),edgept3(2),edgept4(2)]);

width = max([size(im11,2),edgept1(1),edgept2(1),edgept3(1),edgept4(1)]);

t2(height,width,3)=0;

t2 = uint8(t2);

if (ccc~=2)

im11 = curImg2;

end

for i = 1:size(im11,1)

for j = 1:size(im11,2)

t2(i,j,:) = im11(i,j,:);

end

end

p1 = [round(size(im11,2) /2),round(size(im11,1) /2)];

p2 = round(Hfinal*[size(im22,2);size(im22,1);1]);

for i = 1:size(im22,1)

for j = 1:size(im22,2)

x = Hfinal*[j;i;1];

tx = max(1,round(x(2))); ty = max(1,round(x(1)));

if (t2(tx,ty,1)~=0&&t2(tx,ty,2)~=0&&t2(tx,ty,3)~=0)

d1 = sqrt((ty-p1(1))^2+(tx-p1(2))^2);

d2 = sqrt((ty-p2(1))^2+(tx-p2(2))^2);

para1 = d2/(d1+d2);

para2 = d1/(d1+d2);

t2(tx,ty,:) = round((im22(i,j,:).*para2 + t2(tx,ty,:).*para1));

else

t2(tx,ty,:) = im22(i,j,:);

end

end

end

for i =2:height-1

for j = 2:width-1

if (t2(i,j,1)==0&&t2(i,j,2)==0&&t2(i,j,3)==0)

t2(i,j,:) = round(t2(i-1,j-1,:)/ + t2(i-1,j+1,:)/ +t2(i+1,j-1,:)/ +t2(i+1,j+1,:)/ );

end

end

end

%figure;imshow(t2);

curImg2 = t2;

set,’Pos’,[28,28,size(curImg2,2),size(curImg2,1)]);

axes;

imshow(curImg2);