原文來自於http://www.cnblogs.com/naive/p/3579610.html

背景介紹

影象拼接是一項應用廣泛的影象處理技術。根據特徵點的相互匹配，可以將多張小視角的影象拼接成為一張大視角的影象，在廣角照片合成、衛星照片處理、醫學影象處理等領域都有應用。早期的影象拼接主要是運用畫素值匹配的方法。後來，人們分別在兩幅影象中尋找拐點、邊緣等穩定的特徵，用特徵匹配的方法拼接影象。本實驗根據Matthew Brown (2005) 描述的方法，實現多張生活照的拼接。

特徵點捕捉 (Interest Point Detection)

首先，拍攝兩張場景有重合的照片。為了保證有足夠多的公共特徵點，照片的重合度應該保證在30%以上。將兩張照片轉換為灰度影象，對影象做σ=1的高斯模糊。在Matthew的文章中，他建立了一個影象金字塔，在不同尺度尋找Harris關鍵點。考慮到將要拼接的照片視野尺寸接近，故簡化此步驟，僅在原圖提取特徵點。

接下來用sobel運算元計算影象在x、y兩個方向亮度的梯度，用σ=1.5的高斯函式對梯度做平滑處理，減小噪點對亮度的影響。很容易發現，若我們求一小塊區域內亮度的累加值，在影象變化平緩的區域上下左右移動視窗累加值的變化並不明顯；在物體的邊緣，沿著邊緣方向的變化也不明顯；而在關鍵點附近，輕微的移動視窗都會強烈改變亮度的累加值，如圖1所示。

圖1 http://www.cse.psu.edu/~rcollins/CSE486/lecture06.pdf

亮度的變化值可以用下面的公式計算得到：

        (1)

其中，w(x, y) 是高斯函式的權重，I(x, y)是該點亮度的梯度。

在計算時，上面的公式又可以近似為如下：

       (2)

通過比較矩陣的特徵值l1和l2，我們可以判斷該點所處的狀態。若l1>>l2或者l2<<l1，表示該點位於縱向或者橫向的邊緣；若l1和l2近似且值很小，表示該點位於平滑區域；若l1和l2近似但值很大，表示該點位於關鍵點。根據Harris and Stephens (1988) 的介紹，我們並不需要直接計算兩個特徵值，用R = Det(H)/Tr(H)2的值就可以反映兩個特徵值的比值，這樣可以減少運算量。我們保留R > 2的點。除此之外，每個點的R和周圍8鄰域畫素的R值比較，僅保留區域性R值最大的點。最後，去除圖片邊界附近的關鍵點。

至此，我們在兩幅圖片分別得到了一組關鍵點，如圖2所示。

圖2 Harris Corner

自適應非極大值抑制 (Adaptive Non-Maximal Suppression)

由於上一步得到的關鍵點很多，直接計算會導致很大的運算量，也會增加誤差。接下去就要去除其中絕大部分的關鍵點，僅保留一些特徵明顯點，且讓關鍵點在整幅影象內分佈均勻。Matthew發明了adaptive non-maximal suppression (ANMS) 方法來擇優選取特定數量的關鍵點。

ANMS的思想是有一個半徑r，初始值為無限遠。當r不斷減小時，保留在半徑r以內其它關鍵點R值均小於中心點R值的關鍵點，將其加入佇列。佇列內的關鍵點數達到預設值後停止搜尋。

Xi是上一步得到的關鍵點的2維座標，G是所有關鍵點的集合，c=0.9。

實際計算時，我們將上述過程相反。這裡我設定每幅影象各提取500個關鍵點。首先找出整幅圖片R值最大的關鍵點Rmax，加入佇列，並且得到Rmax*0.9的值。遍歷所有關鍵點，若該關鍵點xi的Ri> Rmax*0.9, 該點的半徑設為無限遠；若該關鍵點xi的Ri< Rmax*0.9，計算該點到離它最近的Rj>0.9R的點xi，記錄兩點間的距離ri。最後將所有r排序，找出r最大的500個點，如圖3所示。

圖3 Harris corner after ANMS

關鍵點的描述 (Feature Descriptor)

關鍵點的描述方法有很多種，包括區域性梯度描述、尺度不變特徵變換 (SIFT、SUFT) 等等。因為生活照的旋轉角度通常不超過15°，所以這裡不考慮關鍵點的旋轉不變性。

對影象做適度的高斯模糊，以關鍵點為中心，取40x40畫素的區域。將該區域降取樣至8x8的大小，生成一個64維的向量。對向量做歸一化處理。每個關鍵點都用一個64維的向量表示，於是每幅影象分別得到了一個500x64的特徵矩陣。

關鍵點的匹配

首先，從兩幅圖片的500個特徵點中篩選出配對的點。篩選的方法是先計算500個特徵點兩兩之間的歐氏距離，按照距離由小到大排序。通常情況下選擇距離最小的一對特徵向量配對。Lowe（2004）認為，僅僅觀察最小距離並不能有效篩選配對特徵點，而用最小的距離和第二小的距離的比值可以很好的進行篩選。如圖4所示， 使用距離的比值能夠獲得更高的true positive， 同時控制較低的false positive。我使用的閾值是r1/r2<0.5。經過篩選後的配對特徵點 如圖5所示

圖 4. 配對正確率和配對方法、閾值選擇的關係

圖 5. 篩選後的配對特徵點

關鍵點的匹配使用Random Sample Consensus (RANSAC) 演算法。以一幅影象為基準，每次從中隨機選擇8個點，在另一幅影象中找出配對的8個點。用8對點計算得到一個homography，將基準圖中剩餘的特徵點按照homography變換投影到另一幅影象，統計配對點的個數。

重複上述步驟2000次，得到準確配對最多的一個homography。至此，兩幅影象的投影變換關係已經找到。

新影象的合成

在做影象投影前，要先新建一個空白畫布。比較投影后兩幅影象的2維座標的上下左右邊界，選取各個方向邊界的最大值作為新影象的尺寸。同時，計算得到兩幅影象的交叉區域。

在兩幅影象的交叉區域，按照cross dissolve的方法制作兩塊如圖6所示的蒙版，3個通道的畫素值再次區間內遞減（遞升）。

效果展示

下面展示幾張照片拼接的效果圖。

圖 7. 拼接完成的新影象

圖 8. 以左邊照片為基準拼接

附Matlab程式碼：

function [output_image] = image_stitching(input_A, input_B)
% -------------------------------------------------------------------------
% 1. Load both images, convert to double and to grayscale.
% 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. 
% 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.
% 
% [email protected]

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

% READ IMAGE, GET SIZE INFORMATION
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
[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);

% ADAPTIVE NON-MAXIMAL SUPPRESSION (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);

% 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);
% THE RATIO OF FIRST AND SECOND DISTANCE IS A BETTER CRETIA THAN DIRECTLY
% USING THE DISTANCE. RATIO LESS THAN .5 GIVES AN ACCEPTABLE ERROR RATE.
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);

% USE 4-POINT RANSAC TO COMPUTE A ROBUST HOMOGRAPHY ESTIMATE
% 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,
matcher_B = [y_B, x_B, ones(npoints,1)]'; %!!! so switch x and y here.
[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(:,:,1) = interp2(X, Y, double(image_A(:,:,1)), XX, YY);
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
imshow(uint8(newImage));

% -------------------------------------------------------------------------
% ------------------------------- other functions -------------------------
% -------------------------------------------------------------------------
function [xp, yp, value] = harris(input_image, sigma,thd, r)
% Detect harris corner 
% Input:
% sigma - standard deviation of smoothing Gaussian
% r - radius of region considered in non-maximal suppression
% Output:
% xp - x coordinates of harris corner points
% yp - y coordinates of harris corner points
% value - values of R at harris corner points

% CONVERT RGB IMAGE TO GRAY-SCALE, AND BLUR WITH G1 KERNEL
g1 = fspecial('gaussian', 7, 1);
gray_image = imfilter(input_image, g1);

% FILTER INPUT IMAGE WITH SOBEL KERNEL TO GET GRADIENT ON X AND Y
% ORIENTATION RESPECTIVELY
h = fspecial('sobel');
Ix = imfilter(gray_image,h,'replicate','same');
Iy = imfilter(gray_image,h','replicate','same');

% GENERATE GAUSSIAN FILTER OF SIZE 6*SIGMA (± 3SIGMA) AND OF MINIMUM SIZE 1x1
g = fspecial('gaussian',fix(6*sigma), sigma);

Ix2 = imfilter(Ix.^2, g, 'same').*(sigma^2); 
Iy2 = imfilter(Iy.^2, g, 'same').*(sigma^2);
Ixy = imfilter(Ix.*Iy, g, 'same').*(sigma^2);

% HARRIS CORNER MEASURE
R = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2 + eps); 
% ANOTHER MEASUREMENT, USUALLY k IS BETWEEN 0.04 ~ 0.06
% response = (Ix2.*Iy2 - Ixy.^2) - k*(Ix2 + Iy2).^2;

% GET RID OF CORNERS WHICH IS CLOSE TO BORDER
R([1:20, end-20:end], :) = 0;
R(:,[1:20,end-20:end]) = 0;

% SUPRESS NON-MAX 
d = 2*r+1; 
localmax = ordfilt2(R,d^2,true(d)); 
R = R.*(and(R==localmax, R>thd));

% RETURN X AND Y COORDINATES 
[xp,yp,value] = find(R);

function [newx, newy, newvalue] = ada_nonmax_suppression(xp, yp, value, n)
% Adaptive non-maximun suppression 
% For each Harris Corner point, the minimum suppression radius is the
% minimum distance from that point to a different point with a higher 
% corner strength. 
% Input:
% xp,yp - coordinates of harris corner points
% value - strength of suppression
% n - number of interesting points
% Output:
% newx, newy - new x and y coordinates after adaptive non-maximun suppression
% value - strength of suppression after adaptive non-maximun suppression

% ALLOCATE MEMORY
% newx = zeros(n,1);
% newy = zeros(n,1);
% newvalue = zeros(n,1);

if(length(xp) < n)
newx = xp;
newy = yp;
newvalue = value;
return;
end

radius = zeros(n,1);
c = .9;
maxvalue = max(value)*c;
for i=1:length(xp)
if(value(i)>maxvalue)
radius(i) = 99999999;
continue;
else
dist = (xp-xp(i)).^2 + (yp-yp(i)).^2;
dist((value*c) < value(i)) = [];
radius(i) = sqrt(min(dist));
end
end

[~, index] = sort(radius,'descend');
index = index(1:n);

newx = xp(index);
newy = yp(index);
newvalue = value(index);

function n2 = dist2(x, c)
% DIST2 Calculates squared distance between two sets of points.
% Adapted from Netlab neural network software:
% http://www.ncrg.aston.ac.uk/netlab/index.php
%
% Description
% D = DIST2(X, C) takes two matrices of vectors and calculates the
% squared Euclidean distance between them. Both matrices must be of
% the same column dimension. If X has M rows and N columns, and C has
% L rows and N columns, then the result has M rows and L columns. The
% I, Jth entry is the squared distance from the Ith row of X to the
% Jth row of C.
%
%
% Copyright (c) Ian T Nabney (1996-2001)

[ndata, dimx] = size(x);
[ncentres, dimc] = size(c);
if dimx ~= dimc
error('Data dimension does not match dimension of centres')
end

n2 = (ones(ncentres, 1) * sum((x.^2)', 1))' + ...
ones(ndata, 1) * sum((c.^2)',1) - ...
2.*(x*(c'));

% Rounding errors occasionally cause negative entries in n2
if any(any(n2<0))
n2(n2<0) = 0;
end

function [descriptors] = getFeatureDescriptor(input_image, xp, yp, sigma)
% Extract non-rotation invariant feature descriptors
% Input:
% input_image - input gray-scale image
% xx - x coordinates of potential feature points
% yy - y coordinates of potential feature points
% output:
% descriptors - array of descriptors

% FIRST BLUR WITH GAUSSIAN KERNEL
g = fspecial('gaussian', 5, sigma);
blurred_image = imfilter(input_image, g, 'replicate','same');

% THEN TAKE A 40x40 PIXEL WINDOW AND DOWNSAMPLE TO 8x8 PATCH
npoints = length(xp);
descriptors = zeros(npoints,64);

for i = 1:npoints
%pA = imresize( blurred_image(xp(i)-20:xp(i)+19, yp(i)-20:yp(i)+19), .2);
patch = blurred_image(xp(i)-20:xp(i)+19, yp(i)-20:yp(i)+19);
patch = imresize(patch, .2);
descriptors(i,:) = reshape((patch - mean2(patch))./std2(patch), 1, 64); 
end

function [hh] = getHomographyMatrix(point_ref, point_src, npoints)
% Use corresponding points in both images to recover the parameters of the transformation 
% Input:
% x_ref, x_src --- x coordinates of point correspondences
% y_ref, y_src --- y coordinates of point correspondences
% Output:
% h --- matrix of transformation

% NUMBER OF POINT CORRESPONDENCES
x_ref = point_ref(1,:)';
y_ref = point_ref(2,:)';
x_src = point_src(1,:)';
y_src = point_src(2,:)';

% COEFFICIENTS ON THE RIGHT SIDE OF LINEAR EQUATIONS
A = zeros(npoints*2,8);
A(1:2:end,1:3) = [x_ref, y_ref, ones(npoints,1)];
A(2:2:end,4:6) = [x_ref, y_ref, ones(npoints,1)];
A(1:2:end,7:8) = [-x_ref.*x_src, -y_ref.*x_src];
A(2:2:end,7:8) = [-x_ref.*y_src, -y_ref.*y_src];

% COEFFICIENT ON THE LEFT SIDE OF LINEAR EQUATIONS
B = [x_src, y_src];
B = reshape(B',npoints*2,1);

% SOLVE LINEAR EQUATIONS
h = A\B;

hh = [h(1),h(2),h(3);h(4),h(5),h(6);h(7),h(8),1];

function [hh, inliers] = ransacfithomography(ref_P, dst_P, npoints, threshold);
% 4-point RANSAC fitting
% Input:
% matcher_A - match points from image A, a matrix of 3xN, the third row is 1
% matcher_B - match points from image B, a matrix of 3xN, the third row is 1
% thd - distance threshold
% npoints - number of samples
% 
% 1. Randomly select minimal subset of points
% 2. Hypothesize a model
% 3. Computer error function
% 4. Select points consistent with model
% 5. Repeat hypothesize-and-verify loop
% 
% Yihua Zhao 02-01-2014
% [email protected]

ninlier = 0;
fpoints = 8; %number of fitting points
for i=1:2000
rd = randi([1 npoints],1,fpoints);
pR = ref_P(:,rd);
pD = dst_P(:,rd);
h = getHomographyMatrix(pR,pD,fpoints);
rref_P = h*ref_P;
rref_P(1,:) = rref_P(1,:)./rref_P(3,:);
rref_P(2,:) = rref_P(2,:)./rref_P(3,:);
error = (rref_P(1,:) - dst_P(1,:)).^2 + (rref_P(2,:) - dst_P(2,:)).^2;
n = nnz(error<threshold);
if(n >= npoints*.95)
hh=h;
inliers = find(error<threshold);
pause();
break;
elseif(n>ninlier)
ninlier = n;
hh=h;
inliers = find(error<threshold);
end 
end

function [newH, newW, x1, y1, x2, y2] = getNewSize(transform, h2, w2, h1, w1)
% Calculate the size of new mosaic
% Input:
% transform - homography matrix
% h1 - height of the unwarped image
% w1 - width of the unwarped image
% h2 - height of the warped image
% w2 - height of the warped image
% Output:
% newH - height of the new image
% newW - width of the new image
% x1 - x coordate of lefttop corner of new image
% y1 - y coordate of lefttop corner of new image
% x2 - x coordate of lefttop corner of unwarped image
% y2 - y coordate of lefttop corner of unwarped image
% 
% Yihua Zhao 02-02-2014
% [email protected]
%

% CREATE MESH-GRID FOR THE WARPED IMAGE
[X,Y] = meshgrid(1:w2,1:h2);
AA = ones(3,h2*w2);
AA(1,:) = reshape(X,1,h2*w2);
AA(2,:) = reshape(Y,1,h2*w2);

% DETERMINE THE FOUR CORNER OF NEW IMAGE
newAA = transform\AA;
new_left = fix(min([1,min(newAA(1,:)./newAA(3,:))]));
new_right = fix(max([w1,max(newAA(1,:)./newAA(3,:))]));
new_top = fix(min([1,min(newAA(2,:)./newAA(3,:))]));
new_bottom = fix(max([h1,max(newAA(2,:)./newAA(3,:))]));

newH = new_bottom - new_top + 1;
newW = new_right - new_left + 1;
x1 = new_left;
y1 = new_top;
x2 = 2 - new_left;
y2 = 2 - new_top;

function [newImage] = blend(warped_image, unwarped_image, x, y)
% Blend two image by using cross dissolve
% Input:
% warped_image - original image
% unwarped_image - the other image
% x - x coordinate of the lefttop corner of unwarped image
% y - y coordinate of the lefttop corner of unwarped image
% Output:
% newImage
% 
% Yihua Zhao 02-02-2014
% [email protected]
%

% MAKE MASKS FOR BOTH IMAGES 
warped_image(isnan(warped_image))=0;
maskA = (warped_image(:,:,1)>0 |warped_image(:,:,2)>0 | warped_image(:,:,3)>0);
newImage = zeros(size(warped_image));
newImage(y:y+size(unwarped_image,1)-1, x: x+size(unwarped_image,2)-1,:) = unwarped_image;
mask = (newImage(:,:,1)>0 | newImage(:,:,2)>0 | newImage(:,:,3)>0);
mask = and(maskA, mask);

% GET THE OVERLAID REGION
[~,col] = find(mask);
left = min(col);
right = max(col);
mask = ones(size(mask));
if( x<2)
mask(:,left:right) = repmat(linspace(0,1,right-left+1),size(mask,1),1);
else
mask(:,left:right) = repmat(linspace(1,0,right-left+1),size(mask,1),1);
end

% BLEND EACH CHANNEL
warped_image(:,:,1) = warped_image(:,:,1).*mask;
warped_image(:,:,2) = warped_image(:,:,2).*mask;
warped_image(:,:,3) = warped_image(:,:,3).*mask;

% REVERSE THE ALPHA VALUE
if( x<2)
mask(:,left:right) = repmat(linspace(1,0,right-left+1),size(mask,1),1);
else
mask(:,left:right) = repmat(linspace(0,1,right-left+1),size(mask,1),1);
end
newImage(:,:,1) = newImage(:,:,1).*mask;
newImage(:,:,2) = newImage(:,:,2).*mask;
newImage(:,:,3) = newImage(:,:,3).*mask;

newImage(:,:,1) = warped_image(:,:,1) + newImage(:,:,1);
newImage(:,:,2) = warped_image(:,:,2) + newImage(:,:,2);
newImage(:,:,3) = warped_image(:,:,3) + newImage(:,:,3);