matlab图像拼接1.rar_harris拼接_两图像拼接_角点匹配

% function [output_image] = image_stitching(image_A, image_B) % [output_image] = image_stitching(input_A, input_B) % ------------------------------------------------------------------------- % 1. Load both images, convert to double and to grayscale. % 加载两个图像，转换为double和灰度 % 2. Detect feature points in both images. 在两个图像中检测特征点 % 3. Extract fixed-size patches around every keypoint in both images, and % form descriptors simply by "flattening" the pixel values in each patch to % one-dimensional vectors. % 在两个图像的每个关键点上提取固定大小的补丁，并通过将每个补丁中的像素值“压平”到一维向量来形成描述符 % 4. Compute distances between every descriptor in one image and every descriptor in the other image. % 计算一个图像中的每个描述符与另一个图像中的每个描述符之间的距离 % 5. Select putative matches based on the matrix of pairwise descriptor % distances obtained above. 根据上面得到的成对描述距离矩阵选择假定匹配 % 6. Run RANSAC to estimate (1) an affine transformation and (2) a % homography mapping one image onto the other. % 运行RANSAC来估计(1)仿射变换和(2)一个映射将一个图像映射到另一个 % 7. Warp one image onto the other using the estimated transformation. % 用估计变换将一个图像传送到另一个图像上 % 8. Create a new image big enough to hold the panorama and composite the % two images into it. 创建一个足够大的新图像来容纳全景，并将两个图像合成到它里面 % % Input: % input_A - filename of warped image % input_B - filename of unwarped image % Output: % output_image - combined new image % % Reference: % [1] C.G. Harris and M.J. Stephens, A combined corner and edge detector, 1988. % [2] Matthew Brown, Multi-Image Matching using Multi-Scale Oriented Patches. % % zhyh8341@gmail.com % ------------------------------------------------------------------------- clear clc % READ IMAGE, GET SIZE INFORMATION 读取图像，获取大小信息 input_A=imread('input_A.jpg'); input_B=imread('input_B.jpg'); image_A = input_A; image_B = input_B; % image_A = imread(input_A); % image_B = imread(input_B); [height_wrap, width_wrap,~] = size(image_A); [height_unwrap, width_unwrap,~] = size(image_B); % CONVERT TO GRAY SCALE 转换为灰度 gray_A = im2double(rgb2gray(image_A)); gray_B = im2double(rgb2gray(image_B)); % FIND HARRIS CORNERS IN BOTH IMAGE 在这两个图像中找到harris [x_A, y_A, v_A] = harris(gray_A, 2, 0.0, 2); [x_B, y_B, v_B] = harris(gray_B, 2, 0.0, 2); figure(1) subplot(1,2,1); imshow(gray_A); hold on plot(y_A,x_A,'r+'); subplot(1,2,2); imshow(gray_B); hold on plot(y_B,x_B,'r+'); set(gcf,'Color','w'); % ADAPTIVE NON-MAXIMAL SUPPRESSION (ANMS) 自适应非最大抑制(ANMS) ncorners = 500; [x_A, y_A, ~] = ada_nonmax_suppression(x_A, y_A, v_A, ncorners); [x_B, y_B, ~] = ada_nonmax_suppression(x_B, y_B, v_B, ncorners); figure(2) subplot(1,2,1); imshow(gray_A); hold on plot(y_A,x_A,'r+'); subplot(1,2,2); imshow(gray_B); hold on plot(y_B,x_B,'r+'); set(gcf,'Color','w'); % EXTRACT FEATURE DESCRIPTORS 提取特征描述符 sigma = 7; [des_A] = getFeatureDescriptor(gray_A, x_A, y_A, sigma); [des_B] = getFeatureDescriptor(gray_B, x_B, y_B, sigma); % IMPLEMENT FEATURE MATCHING 实现特征匹配 dist = dist2(des_A,des_B); [ord_dist, index] = sort(dist, 2); % 自带排序函数 % [Y,I] = sort(X,DIM,MODE) % sort函数默认Mode为'ascend'为升序，sort(X,'descend')为降序排列。 % sort(X)若X是矩阵，默认对X的各列进行升序排列 % sort(X,dim) % dim=1时等效sort(X) % dim=2时表示对X中的各行元素升序排列 % 若欲保留排列前的索引，则可用[sX,index] = sort(X) ，排序后，sX是排序好的向量，index是 向量sX中对X 的索引 % THE RATIO OF FIRST AND SECOND DISTANCE IS A BETTER CRETIA THAN DIRECTLY % USING THE DISTANCE. 第一和第二距离的比值比直接使用距离更好 % RATIO LESS THAN 0.5 GIVES AN ACCEPTABLE ERROR RATE. 小于0.5的比率给出一个可接受的错误率 ratio = ord_dist(:,1)./ord_dist(:,2); threshold = 0.5; idx = ratio<threshold; x_A = x_A(idx); y_A = y_A(idx); x_B = x_B(index(idx,1)); y_B = y_B(index(idx,1)); npoints = length(x_A); figure(3) subplot(1,2,1); imshow(gray_A); hold on plot(y_A,x_A,'r+'); subplot(1,2,2); imshow(gray_B); hold on plot(y_B,x_B,'r+'); set(gcf,'Color','w'); % USE 4-POINT RANSAC TO COMPUTE A ROBUST HOMOGRAPHY ESTIMATE % 使用4点RANSAC来计算一个强健的同调估计 % KEEP THE FIRST IMAGE UNWARPED, WARP THE SECOND TO THE FIRST % 保持第一个图像不变形，将第二个图像扭曲到第一个 matcher_A = [y_A, x_A, ones(npoints,1)]'; %!!! previous x is y and y is x, 之前的x是y, y是x matcher_B = [y_B, x_B, ones(npoints,1)]'; %!!! so switch x and y here. 把x和y交换一下 [hh, ~] = ransacfithomography(matcher_B, matcher_A, npoints, 10); % s = load('matcher.mat'); % matcher_A = s.matcher(1:3,:); % matcher_B = s.matcher(4:6,:); % npoints = 60; % [hh, inliers] = ransacfithomography(matcher_B, matcher_A, npoints, 10); % USE INVERSE WARP METHOD 使用逆经法 % DETERMINE THE SIZE OF THE WHOLE IMAGE 确定整个图像的大小 [newH, newW, newX, newY, xB, yB] = getNewSize(hh, height_wrap, width_wrap, height_unwrap, width_unwrap); [X,Y] = meshgrid(1:width_wrap,1:height_wrap); % 生成二、三维网格 [XX,YY] = meshgrid(newX:newX+newW-1, newY:newY+newH-1); AA = ones(3,newH*newW); AA(1,:) = reshape(XX,1,newH*newW); % 重新按列排列 AA(2,:) = reshape(YY,1,newH*newW); AA = hh*AA; XX = reshape(AA(1,:)./AA(3,:), newH, newW); YY = reshape(AA(2,:)./AA(3,:), newH, newW); % INTERPOLATION, WARP IMAGE A INTO NEW IMAGE 插值，使图像成为新的图像 % newImage=[:,:,3]; newImage(:,:,1) = interp2(X, Y, double(image_A(:,:,1)), XX, YY); % 插值函数 % ZI = interp2(X,Y,Z,XI,YI) % X,Y是原始数据，相当于坐标，类似于meshgrid的坐标范围，这么说应该很容易理解...... % Z是在上述坐标下的数值，也就是在坐标[xi yi]下的zi % XI,YI就是用于插值的坐标，返回值ZI就是用于提取插值之后，对应位置的值 % 这里需要注意： % X 与Y必须是单调的若Xi与Yi中有在X与Y范围之外的点，则相应地返回nan(Not a Number) % ZI = interp2(Z,XI,YI) % 缺省地，X=1:n、Y=1:m，其中[m,n]=size(Z)。再按第一种情形进行计算 % ZI = interp2(Z,n) % 作n次递归计算，在Z的每两个元素之间插入它们的二维插值，这样，Z的阶数将不断增加。 % interp2(Z)等价于interp2(z,1) % ZI = interp2(X,Y,Z,XI,YI,method) % 用指定的算法method 计算二维插值： % ’linear’：双线性插值算法(缺省算法); % ’nearest’：最临近插值; % ’spline’：三次样条插值; % ’cubic’：双三次插值 newImage(:,:,2) = interp2(X, Y, double(image_A(:,:,2)), XX, YY); newImage(:,:,3) = interp2(X, Y, double(image_A(:,:,3)), XX, YY); % BLEND IMAGE BY CROSS DISSOLVE 混合图像通过交叉溶解 [newImage] = blend(newImage, image_B, xB, yB); % DISPLAY IMAGE MOSIAC 显示图象MOSIAC figure(5); imshow(uint8(newImage));