用matlab编写一个多幅图片融合拼接程序_百度知道
用matlab编写一个多幅图片融合拼接程序
用matlab编写一个多幅图片融合拼接程序有偿...
用matlab编写一个多幅图片融合拼接程序有偿
答题抽奖
首次认真答题后
即可获得3次抽奖机会，100%中奖。
liu_wei128
liu_wei128
采纳数：98
获赞数：75
什么意思。把几张图片拼在一起吗
为你推荐：
其他类似问题
个人、企业类
违法有害信息,请在下方选择后提交
色情、暴力
我们会通过消息、邮箱等方式尽快将举报结果通知您。本博客与以下文档资料一起服用效果更佳。
Matlab源码地址：
开始正文。
梳理一下本篇博客图像拼接的原理：
特征检测：SIFT角点检测
特征描述：SIFT描述子
特征匹配：暴力搜索+欧氏距离
求变换矩阵：最小二乘法+RANSAC+仿射变换
拼接多幅图像
完整实现一个SIFT都很麻烦，也并没必要。本博客是在斯坦福大学计算机视觉课程大作业的基础上实现的，读者也可以按课程PPT及作业pdf指导书，自己实现一遍，定会加深自己对其中关键算法的理解。
本博客根据其作业指导书，展示其中一些算法的Matlab代码实现。
全景拼接是计算机视觉领域取得的一项早期成就。在2007年，Brown 和 Lowe发表了著名的图像拼接论文。自打那以后，自动全景拼接技术受到广泛应用，例如Google街景地图、智能手机上的全景照片、拼接软件比如Photosynth和AutoStitch。
在这个大作业中，我们会从多幅图像中匹配SIFT关键点，来构建单幅全景图片。这涉及以下任务：
使用高斯差分（DoG）检测器找关键点，返回它的坐标位置和尺度。（这一步已经提供给你）
给一副图像的每个关键点构建SIFT描述子。
从两幅不同的图像中比较两组SIFT描述子，找到匹配点。
给定一个关键点匹配的列表，使用最小二乘法找到仿射变换矩阵，这个矩阵能将iamge1上的位置映射到image2上的位置。
使用RANSAC使仿射变换矩阵的估计具有更好的鲁棒性。
给定变换矩阵，使用它变换（平移、尺度、或者倾斜）image1，将它覆盖到image2上面，构建一个全景图。（这一步已经提供给你）
在真实世界场景中的特定例子中，把多幅图像拼接在一起。
2.构建SIFT描述子
复习SIFT算法的讲义PPT，编写给定的SIFTDescriptor.m，来为每个DoG关键点，产生SIFT关键点描述子。注意关键点的位置和尺度已经提供给你，使用的是高斯金字塔。
运行提供的EvaluateSIFTDescriptor.m，检查你的实现。
这一步需要自己写的代码包括计算图像导数（x和y方向）、梯度（大小和方向）、计算邻域块的主方向、每个采样点的梯度相对方向、计算直方图的变量（方向箱子的划分，每个箱子的幅值）、计算梯度直方图以及将其串接进一个1*128的数组中。
代码提供了标准的结果数据和测试代码，经过测试，我编写的SIFT与标准数据偏差明显。
注：感谢yuanyuan12222的评论指正。出现偏差的原因是：主方向选取没有向下取（308行注释有说明），将312行代码direction =2*pi*loc(1)/num_改为direction =2*pi*(loc(1)-1)/num_就OK了。
说明我的某部分实现不符合标准实现。但由最后的拼接效果来看，实现是可用的。
我的实现（Matlab代码）：
SIFTDescriptor.m
function descriptors = SIFTDescriptor(pyramid, keyPtLoc, keyPtScale)
% SIFTDescriptor Build SIFT descriptors from image at detected key points'
% location with detected key points' scale and angle
pyramid: Image pyramid. pyramid{i} is a rescaled version of the
original image, in grayscale double format
keyPtLoc: N * 2 matrix, each row is a key point pixel location in
pyramid{round(keyPtScale)}. So pyramid{round(keyPtScale)}(y,x) is the center of the keypoint
keyPtScale: N * 1 matrix, each entry holds the index in the Gaussian
pyramid for the keypoint. Earlier code interpolates things, so this is a
double, but we'll just round it to an integer.
descriptors: N * 128 matrix, each row is a feature descriptor
% Precompute the gradients at all pixels of all pyramid scales
% This is a cell array, which is like an ArrayList that holds matrices.
% You use {} to index into it, like this: magnitude = grad_mag{1}
grad_theta = cell(length(pyramid),1);
grad_mag = cell(length(pyramid),1);
for scale = 1:length(pyramid)
currentImage = pyramid{scale};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
YOUR CODE HERE
Read the doc for filter2.
Use with the filter [-1 0 1] to fill in img_dx,
and the filter [-1;0;1] to fill in img_dy.
Please use the filter2 'same' option so
the result will be the same size as the image.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% gradient image, for gradients in x direction.
img_dx = zeros(size(currentImage));
img_dx=filter2([-1 0 1],currentImage);
% gradients in y direction.
img_dy = zeros(size(currentImage));
img_dy=filter2([-1;0;1],currentImage);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
END OF YOUR CODE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
YOUR CODE HERE
Use img_dx and img_dy to compute the magnitude
and angle of the gradient at each pixel.
store them in grad_mag{scale} and grad_theta{scale} respectively.
The atan2 function will be helpful for calculating angle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the magnitude and orientation of the gradient.
grad_mag{scale} = zeros(size(currentImage));
grad_theta{scale} = zeros(size(currentImage));
grad_mag{scale}=sqrt(img_dx.^2+img_dy.^2);
grad_theta{scale}=atan2(img_dy,img_dx);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
END OF YOUR CODE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% atan2 gives angles from -pi to pi. To make the histogram code
% easier, we'll change that to 0 to 2*pi.
grad_theta{scale} = mod(grad_theta{scale}, 2*pi);
% The number of bins into which each gradient vector will be placed.
num_angles = 8;
% The patch extracted around each keypoint will be divided into a grid
% of num_histograms x num_histograms.
num_histograms = 4;
% Each histogram covers an area "pixelsPerHistogram" wide and tall
pixelsPerHistogram = 4;
% For each keypoint we will extract a region of size
% patch_size x patch_size centered at the keypoint.
patch_size = num_histograms * pixelsPerH
% Number of keypoints that were found by the DoG blob detector
N = size(keyPtLoc, 1);
% Initialize descriptors to zero
descriptors = zeros(N, num_histograms * num_histograms * num_angles);
% Iterate over all keypoints
for i = 1 : N
scale = round(keyPtScale(i));
% Find the window of pixels that contributes to the descriptor for the
% current keypoint.
xAtScale = keyPtLoc(i, 1);%center of the DoG keypoint in the pyramid{2} image
yAtScale = keyPtLoc(i, 2);
x_lo = round(xAtScale - patch_size / 2);
x_hi = x_lo+patch_size-1;
y_lo = round(yAtScale - patch_size / 2);
y_hi = y_lo+patch_size-1;
% These are the gradient magnitude and angle images from the
% correct scale level. You computed these above.
magnitudes = grad_mag{scale};
thetas = grad_theta{scale};
% Extract the patch from that window around the keypoint
patch_mag = zeros(patch_size,patch_size);
patch_theta = zeros(patch_size,patch_size);
patch_mag = magnitudes(y_lo:y_hi,x_lo:x_hi);
patch_theta = thetas(y_lo:y_hi,x_lo:x_hi);
% If any keypoint is too close to the boundary of the image
% then we just skip it.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
YOUR CODE HERE:
Express gradient directions relative to the dominant gradient direction
of this keypoint.
HINT: Use the ComputeDominantDirection function below.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 1: compute the dominant gradient direction of the patch
patch_angle_offset = ComputeDominantDirection(patch_mag, patch_theta);
% Step 2: change patch_theta so it's relative to the dominant direction
patch_theta = patch_theta-patch_angle_
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
END OF YOUR CODE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This line will re-map patch_theta into the range 0 to 2*pi
patch_theta = mod(patch_theta, 2*pi);
% Weight the gradient magnitudes using a gaussian function
patch_mag = patch_mag .* fspecial('gaussian', patch_size, patch_size / 2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
YOUR CODE HERE:
Compute the gradient histograms and concatenate them in the
feature variable to form a size 1x128 SIFT descriptor for this keypoint.
HINT: Use the ComputeGradientHistogram function below.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The patch we have extracted should be subdivided into
% num_histograms x num_histograms cells, each of which is size
% pixelsPerHistogram x pixelsPerHistogram.
% Compute a gradient histogram for each cell, and concatenate
% the histograms into a single feature vector of length 128.
% Please traverse the patch row by row, starting in the top left,
% in order to match the given solution. E.g. if you use two
% nested 'for' loops, the loop over x should
be the inner loop.
% (Note: Unlike the SIFT paper, we will not smooth a gradient across
% nearby histograms. For simplicity, we will just assign all
% gradient pixels within a pixelsPerHistogram x pixelsPerHistogram
% square to the same histogram.)
% Initializing the feature vector to size 0. Hint: you can
% concatenate the histogram descriptors to it like this:
% feature = [feature, histogram]
feature = [];
subdivided_patch_theta=zeros(pixelsPerHistogram,pixelsPerHistogram);
subdivided_patch_mag=zeros(pixelsPerHistogram,pixelsPerHistogram);
for y=1:4:13
for x=1:4:13
subdivided_patch_theta=patch_theta(y:y+3,x:x+3);
subdivided_patch_mag=patch_mag(y:y+3,x:x+3);
[histogram,angles]=ComputeGradientHistogram(num_angles,subdivided_patch_mag,subdivided_patch_theta);
feature=[feature,histogram];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
END YOUR CODE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Add the feature vector we just computed to our matrix of SIFT
% descriptors.
descriptors(i, :) =
% Normalize the descriptors.
descriptors = NormalizeDescriptors(descriptors);
function [histogram, angles] = ComputeGradientHistogram(num_bins, gradient_magnitudes, gradient_angles)
% Compute a gradient histogram using gradient magnitudes and directions.
% Each point is assigned to one of num_bins depending on its gradient
% the gradient magnitude of that point is added to its bin.
% num_bins: The number of bins to which points should be assigned.
% gradient_magnitudes, gradient angles:
Two arrays of the same shape where gradient_magnitudes(i) and
gradient_angles(i) give the magnitude and direction of the gradient
for the ith point. gradient_angles ranges from 0 to 2*pi
% histogram: A 1 x num_bins array containing the gradient histogram. Entry 1 is
the sum of entries in gradient_magnitudes whose corresponding
gradient_angles lie between 0 and angle_step. Similarly, entry 2 is for
angles between angle_step and 2*angle_step. Angle_step is calculated as
2*pi/num_bins.
% angles: A 1 x num_bins array which holds the histogram bin lower bounds.
In other words, histogram(i) contains the sum of the
gradient magnitudes of all points whose gradient directions fall
in the range [angles(i), angles(i + 1))
angle_step = 2 * pi / num_
angles = 0 : angle_step : (2*pi-angle_step);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
YOUR CODE HERE:
Use the function inputs to calculate the histogram variable,
as defined above.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
histogram = zeros(1, num_bins);
[rows,cols]=size(gradient_magnitudes);
for m=1:rows
for n=1:cols
for angle=0:angle_step:(2*pi-angle_step)
if(angle&=gradient_angles(m,n)&&gradient_angles(m,n)&angle+angle_step)
histogram(round(angle/angle_step)+1)=histogram(round(angle/angle_step)+1)+gradient_magnitudes(m,n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
END OF YOUR CODE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function direction = ComputeDominantDirection(gradient_magnitudes, gradient_angles)
% Computes the dominant gradient direction for the region around a keypoint
% given the scale of the keypoint and the gradient magnitudes and gradient
% angles of the pixels in the region surrounding the keypoint.
% gradient_magnitudes, gradient_angles:
Two arrays of the same shape where gradient_magnitudes(i) and
gradient_angles(i) give the magnitude and direction of the gradient for
the ith point.
% scale: The scale of the keypoint
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
YOUR CODE HERE:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute a gradient histogram using the weighted gradient magnitudes.
% In David Lowe's paper he suggests using 36 bins for this histogram.
num_bins = 36;
% compute the 36-bin histogram of angles using ComputeGradientHistogram()
[histogram, angles] = ComputeGradientHistogram(num_bins, gradient_magnitudes, gradient_angles);
% Find the maximum value of the gradient histogram, and set "direction"
% to the angle corresponding to the maximum. (To match our solutions,
% just use the lower-bound angle of the max histogram bin. (E.g. return
% 0 radians if it's bin 1.)
peak=max(histogram);
loc=find(histogram==peak);
direction =2*pi*loc(1)/num_
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
END YOUR CODE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function descriptors = NormalizeDescriptors(descriptors)
% Normalizes SIFT descriptors so they're unit vectors. You don't need to
% edit this function.
% descriptors: N x 128 matrix where each row is a SIFT descriptor.
% descriptors: N x 128 matrix containing a normalized version of the input.
% normalize all descriptors so they become unit vectors
lengths = sqrt(sum(descriptors.^2, 2));
nonZeroIndices = find(lengths);
lengths(lengths == 0) = 1;
descriptors = descriptors ./ repmat(lengths, [1 size(descriptors,2)]);
% suppress large entries
descriptors(descriptors & 0.2) = 0.2;
% finally, renormalize to unit length
lengths = sqrt(sum(descriptors.^2, 2));
lengths(lengths == 0) = 1;
descriptors = descriptors ./ repmat(lengths, [1 size(descriptors,2)]);
3.匹配SIFT描述子
编写SIFTSimpleMatcher.m，给定一个image1中SIFT描述子和所有image2中的SIFT描述子，计算它们之间的欧式距离。然后使用它来确定是否是良好的匹配：如果与最近邻的向量的距离，比次近邻向量小的多，我们就认为它是一对匹配。输出结果是一个数组，每行两个值，分别代表匹配的描述子的索引编号。
运行提供的EvaluateSIFTMatcher.m来检查你的实现。
这一步的任务比较单纯，计算两个向量之间的欧式距离，将结果向量排序等。测试结果：
SIFTSimpleMatcher.m
function match = SIFTSimpleMatcher(descriptor1, descriptor2, thresh)
% SIFTSimpleMatcher
Match one set of SIFT descriptors (descriptor1) to another set of
descriptors (decriptor2). Each descriptor from descriptor1 can at
most be matched to one member of descriptor2, but descriptors from
descriptor2 can be matched more than once.
Matches are determined as follows:
For each descriptor vector in descriptor1, find the Euclidean distance
between it and each descriptor vector in descriptor2. If the smallest
distance is less than thresh*(the next smallest distance), we say that
the two vectors are a match, and we add the row [d1 index, d2 index] to
the "match" array.
descriptor1: N1 * 128 matrix, each row is a SIFT descriptor.
descriptor2: N2 * 128 matrix, each row is a SIFT descriptor.
thresh: a given threshold of ratio. Typically 0.7
Match: N * 2 matrix, each row is a match.
For example, Match(k, :) = [i, j] means i-th descriptor in
descriptor1 is matched to j-th descriptor in descriptor2.
if ~exist('thresh', 'var'),
thresh = 0.7;
match = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
YOUR CODE HERE:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[N1,~]=size(descriptor1);
[N2,~]=size(descriptor2);
for i=1:N1
distance=[];
for j=1:N2
subtract=descriptor1(i,:)-descriptor2(j,:);
distance=[distance,norm(subtract)];
sort_distance=sort(distance);
if(sort_distance(1)&0.7*sort_distance(2))
j=find(distance==sort_distance(1));
match=[[i,j]];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
END YOUR CODE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4.计算变换矩阵
我们现在有一个两幅图像之间匹配点的列表。我们会使用它来找到一个变换矩阵，这个矩阵能将image1的点映射到image2对应的坐标系。换句话说，如果在image1中的点x1y1匹配image2中的点x2y2，我们需要找到一个变换矩阵H
???x2y21???=H???x1y11???
需要有足够多的点，Matlab能够帮我们计算最佳的H。复习以前的有关最小二乘法的作业。编写ComputeAffineMatrix.m，根据给定的匹配点对计算H。
运行提供的EvaluateAffineMatrix.m来检查你的实现。
采用最小二乘法求取矩阵方程的解，在Matlab上一个“反斜杠”就能够解决。
ComputeAffineMatrix.m
function H = ComputeAffineMatrix( Pt1, Pt2 )
%ComputeAffineMatrix
Computes the transformation matrix that transforms a point from
coordinate frame 1 to coordinate frame 2
Pt1: N * 2 matrix, each row is a point in image 1
(N must be at least 3)
Pt2: N * 2 matrix, each row is the point in image 2 that
matches the same point in image 1 (N should be more than 3)
H: 3 * 3 affine transformation matrix,
such that H*pt1(i,:) = pt2(i,:)
N = size(Pt1,1);
if size(Pt1, 1) ~= size(Pt2, 1),
error('Dimensions unmatched.');
elseif N&3
error('At least 3 points are required.');
% Convert the input points to homogeneous coordintes.
P1 = [Pt1';ones(1,N)];
P2 = [Pt2';ones(1,N)];
% Now, we must solve for the unknown H that satisfies H*P1=P2
% But MATLAB needs a system in the form Ax=b, and A\b solves for x.
% In other words, the unknown matrix must be on the right.
% But we can use the properties of matrix transpose to get something
% in that form. Just take the transpose of both sides of our equation
% above, to yield P1'*H'=P2'. Then MATLAB can solve for H', and we can
% transpose the result to produce H.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
YOUR CODE HERE:
Use MATLAB's "A\b" syntax to solve for H_transpose as discussed
above, then convert it to the final H
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H_transpose=P1'\P2';
H=H_transpose';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
END OF YOUR CODE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sometimes numerical issues cause least-squares to produce a bottom
% row which is not exactly [0 0 1], which confuses some of the later
% code. So we'll ensure the bottom row is exactly [0 0 1].
H(3,:) = [0 0 1];
我们使用RANSAC（“RANdom SAmple Consensus”）选取内点来计算单应矩阵，而不是直接把所有的SIFT关键点匹配放进ComputeAffineMatrix.m产生结果。在这里，内点是指符合相同的变换矩阵的匹配对。我们已经替你实现RANSAC，除了判断两个点适应矩阵H的符合度的函数没有给。在RANSAC.m里编写ComputeError()函数，计算Hp1和p2之间的欧式距离。
其中，∥∥2代表欧式距离，正如上面第三部分定义的一样。在你完成RANSANCFit.m之后，你可以运行TransformationTester.m来测试你的代码。你应该得到跟图1类似的结果。
我的实现：
RANSACFit.m
function H = RANSACFit(p1, p2, match, maxIter, seedSetSize, maxInlierError, goodFitThresh )
%RANSACFit Use RANSAC to find a robust affine transformation
p1: N1 * 2 matrix, each row is a point
p2: N2 * 2 matrix, each row is a point
match: M * 2 matrix, each row represents a match [index of p1, index of p2]
maxIter: the number of iterations RANSAC will run
seedNum: The number of randomly-chosen seed points that we'll use to fit
our initial circle
maxInlierError: A match not in the seed set is considered an inlier if
its error is less than maxInlierError. Error is
measured as sum of Euclidean distance between transformed
point1 and point2. You need to implement the
ComputeCost function.
goodFitThresh: The threshold for deciding whether or not a model is
for a model to be good, at least goodFitThresh
non-seed points must be declared inliers.
H: a robust estimation of affine transformation from p1 to p2
N = size(match, 1);
error('not enough matches to produce a transformation matrix')
if ~exist('maxIter', 'var'),
maxIter = 200;
if ~exist('seedSetSize', 'var'),
seedSetSize = ceil(0.2 * N);
seedSetSize = max(seedSetSize,3);
if ~exist('maxInlierError', 'var'),
maxInlierError = 30;
if ~exist('goodFitThresh', 'var'),
goodFitThresh = floor(0.7 * N);
H = eye(3);
% below is an obfuscated version of RANSAC. You don't need to
% edit any of this code, just the ComputeError() function below
iota = Inf;
kappa = 0;
alpha = seedSetS
for i = 1 : maxIter,
[beta, gamma] = part(match, alpha);
eta = ComputeAffineMatrix(p1(beta(:, 1), :), p2(beta(:, 2), :));
delta = ComputeError(eta, p1, p2, gamma);
epsilon = (delta &= maxInlierError);
if sum(epsilon(:)) + alpha &= goodFitThresh,
zeta = [ gamma(epsilon, :)];
eta = ComputeAffineMatrix(p1(zeta(:, 1), :), p2(zeta(:, 2), :));
theta = sum(ComputeError(eta, p1, p2, zeta));
if theta & iota,
if sum(sum((H - eye(3)).^2)) == 0,
disp('No RANSAC fit was found.')
function dists = ComputeError(H, pt1, pt2, match)
% Compute the error using transformation matrix H to
% transform the point in pt1 to its matching point in pt2.
H: 3 x 3 transformation matrix where H * [x; 1] transforms the point
(x, y) from the coordinate system of pt1 to the coordinate system of
pt1: N1 x 2 matrix where each ROW is a data point [x_i, y_i]
pt2: N2 x 2 matrix where each ROW is a data point [x_i, y_i]
match: M x 2 matrix, each row represents a match [index of pt1, index of pt2]
dists: An M x 1 vector where dists(i) is the error of fitting the i-th
match to the given transformation matrix.
Error is measured as the Euclidean distance between (transformed pt1)
and pt2 in homogeneous coordinates.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
YOUR CODE HERE.
Convert the points to a usable format, perform the
transformation on pt1 points, and find their distance to their
MATCHING pt2 points.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% hint: If you have an array of indices, MATLAB can directly use it to
% index into another array. For example, pt1(match(:, 1),:) returns a
% matrix whose first row is pt1(match(1,1),:), second row is
% pt1(match(2,1),:), etc. (You may use 'for' loops if this is too
% confusing, but understanding it will make your code simple and fast.)
dists = zeros(size(match,1),1);
transform_pt1=H*[pt1(match(:,1),:)';ones(1,size(match,1))];
subtract=pt2(match(:,2),:)-transform_pt1(1:2,:)';
dists=sqrt(subtract(:,1).^2+subtract(:,2).^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
END YOUR CODE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if size(dists,1) ~= size(match,1) || size(dists,2) ~= 1
error('wrong format');
function [D1, D2] = part(D, splitSize)
idx = randperm(size(D, 1));
D1 = D(idx(1:splitSize), :);
D2 = D(idx(splitSize+1:end), :);
6.拼接多幅图像
我们提供了一个函数，使用你之前写的代码就能高效的拼接一个有序的图像序列（多幅图像是在一条线上拍摄的，每张图像都是最终全景的一部分）。给定一个包含m张图像的序列，
Img1，Img2，……Imgm，
我们的代码使用每相邻两幅图像，然后计算它们之间的变换矩阵，比如讲一个点从Imgi的坐标系转换成Imgi+1的坐标系。（是通过在每对图像上简单的调用你写的代码实现它的。）
我们之后选取一个参考图像Imgref，它处在矩阵序列的中间。我们想让我们最终的全景图处在Imgref的坐标系下。所以，对于对于非参考图像的Imgi，我们需要一个变换矩阵，来将第i帧图像的点转换到ref帧上。（MATLAB然后能够使用这个变换矩阵，帮我们变换图像。）
你的任务是在MultipleStitch.m中，实现函数makeTransformToReferenceFrame。给你一个矩阵列表，包含i帧到i+1帧的转换关系。你必须使用这些矩阵产生一个转换给定帧到参考帧的变换矩阵。
完成了这个部分，你可以通过运行StitchTester.m来检查你的代码。你会得到跟图2相似的结果……
这个任务也相对单纯，只是用到矩阵的乘法和求逆运算。
MultipleStitch.m
function Pano = MultipleStitch( IMAGES, TRANS, fileName )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%MultipleStitch
This function stitches multiple images together and outputs the panoramic stitched image
with a chain of input images and its corresponding transformations.
Given a chain of images:
I1 -& I2 -& I3 -& ... -& Im
and its corresponding transformations:
T1 transforms I1 to I2
T2 transforms I2 to I3
Tm-1 transforms Im-1 to Im
We choose the middle image as the reference image, and the outputed
panorama is in the same coordinate system as the reference image.
For this part, all the image stitching code has been provided to you.
The main task for you is to fill in the code so that current
transformations are used when we produce the final panorama.
Originally, we have
I1 -& I2 -& ... -& Iref -& ... -& Im-1 -& Im
When we fix Iref as the final coordinate system, we want all other
images transformed to Iref. You are responsible for finding the current
transformations used under this circumstances.
IMAGES: 1 * m cell array, each cell contains an image
TRANS: 1 * (m-1) cell array, each cell i contains an affine
transformation matrix that transforms Ii to Ii+1.
fileName: the output file name.
Pano: the final panoramic image.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~exist('fileName', 'var'),
fileName = 'pano.jpg';
if length(IMAGES) ~= length(TRANS)+1,
error('Number of images does not match the number of transformations.');
%% Outbounds of panorama image
outBounds = zeros(2,2);
outBounds(1,:) = Inf;
outBounds(2,:) = -Inf;
%% Choose reference image Iref
refIdx = ceil(median(1:length(IMAGES)));
%% Estimate the largest possible panorama size
[nrows, ncols, ~] = size(IMAGES{1});
nrows = length(IMAGES) *
ncols = length(IMAGES) *
% imageToRefTrans is a 1 x m cell array where imageToRefTrans{i} gives the
% affine transformation from IMAGES{i} to the reference image
% IMAGES{refIdx}. Your task is to fill in this array.
imageToRefTrans = cell(1, length(IMAGES));
% Initialize imageToRefTrans to contain the identity transform.
for idx = 1:length(imageToRefTrans)
imageToRefTrans{idx} = eye(3);
%% Find the correct transformations used for images on the left side of Iref
for idx = refIdx-1:-1:1,
imageToRefTrans{idx} = makeTransformToReferenceFrame(TRANS, idx, refIdx);
T = imageToRefTrans{idx};
tmpBounds = findbounds(maketform('affine', T'), [1 1; ncols nrows]);
outBounds(1,:) = min(outBounds(1,:),tmpBounds(1,:));
outBounds(2,:) = max(outBounds(2,:),tmpBounds(2,:));
%% Find the correct transformations used for images on the right side of Iref
for idx = refIdx + 1 : length(imageToRefTrans),
imageToRefTrans{idx} = makeTransformToReferenceFrame(TRANS, idx, refIdx);
T = imageToRefTrans{idx};
T(3, :) = [0, 0, 1]; % Fix rounding errors in the last row.
tmpBounds = findbounds(maketform('affine', T'), [1 1; ncols nrows]);
outBounds(1,:) = min(outBounds(1,:),tmpBounds(1,:));
outBounds(2,:) = max(outBounds(2,:),tmpBounds(2,:));
%% Stitch the Iref image.
XdataLimit = round(outBounds(:,1)');
YdataLimit = round(outBounds(:,2)');
Pano = imtransform( im2double(IMAGES{refIdx}), maketform('affine', eye(3)), 'bilinear', ...
'XData', XdataLimit, 'YData', YdataLimit, ...
'FillValues', NaN, 'XYScale',1);
%% Transform the images from the left side of Iref using the correct transformations you computed
for idx = refIdx-1:-1:1,
T = imageToRefTrans{idx};
Tform = maketform('affine', T');
AddOn = imtransform(im2double(IMAGES{idx}), Tform, 'bilinear', ...
'XData', XdataLimit, 'YData', YdataLimit, ...
'FillValues', NaN, 'XYScale',1);
result_mask = ~isnan(Pano(:,:,1));
temp_mask = ~isnan(AddOn(:,:,1));
add_mask = temp_mask & (~result_mask);
for c = 1 : size(Pano,3),
cur_im = Pano(:,:,c);
temp_im = AddOn(:,:,c);
cur_im(add_mask) = temp_im(add_mask);
Pano(:,:,c) = cur_
%% Transform the images from the right side of Iref using the correct transformations you computed
for idx = refIdx + 1 : length(imageToRefTrans),
T = imageToRefTrans{idx};
T(3, :) = [0, 0, 1]; % Fix rounding errors in the last row.
Tform = maketform('affine', T');
AddOn = imtransform(im2double(IMAGES{idx}), Tform, 'bilinear', ...
'XData', XdataLimit, 'YData', YdataLimit, ...
'FillValues', NaN, 'XYScale',1);
result_mask = ~isnan(Pano(:,:,1));
temp_mask = ~isnan(AddOn(:,:,1));
add_mask = temp_mask & (~result_mask);
for c = 1 : size(Pano,3),
cur_im = Pano(:,:,c);
temp_im = AddOn(:,:,c);
cur_im(add_mask) = temp_im(add_mask);
Pano(:,:,c) = cur_
%% Cropping the final panorama to leave out black spaces.
[I, J] = ind2sub([size(Pano, 1), size(Pano, 2)], find(~isnan(Pano(:, :, 1))));
upper = max(min(I)-1, 1);
lower = min(max(I)+1, size(Pano, 1));
left = max(min(J)-1, 1);
right = min(max(J)+1, size(Pano, 2));
Pano = Pano(upper:lower, left:right,:);
imwrite(Pano, fileName);
function T = makeTransformToReferenceFrame(i_To_iPlusOne_Transform, currentFrameIndex, refFrameIndex)
%makeTransformToReferenceFrame
i_To_iPlusOne_Transform: this is a cell array where
i_To_iPlusOne_Transform{i} contains the 3x3 homogeneous transformation
matrix that transforms a point in frame i to the corresponding point in
currentFrameIndex: index of the current coordinate frame in i_To_iPlusOne_Transform
refFrameIndex: index of the reference coordinate frame
T: A 3x3 homogeneous transformation matrix that would convert a point
in the current frame into the corresponding point in the reference
frame. For example, if the current frame is 2 and the reference frame
is 3, then T = i_To_iPlusOne_Transform{2}. If the current frame and
reference frame are not adjacent, T will need to be calculated.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
YOUR CODE HERE: Calculate T as defined above.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HINT 1: There are two separate cases to consider: currentFrameIndex &
% refFrameIndex (this is the easier case), and currentFrameIndex &
% refFrameIndex (this is the harder case).
% HINT 2: You can use the pinv function to invert a transformation.
if currentFrameIndex&refFrameIndex
for i=currentFrameIndex:refFrameIndex-1
T=T*i_To_iPlusOne_Transform{i};
elseif currentFrameIndex&refFrameIndex
for i=currentFrameIndex-1:refFrameIndex
T=T*pinv(i_To_iPlusOne_Transform{i});
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
END OF YOUR CODE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
最终测试：
拼接全景结果
全景视频拼接（四）：循环将两幅图像拼接为全景图片
项目要求：利用双摄像头同时采集两个视频，离线拼接，将两个视频拼接成一个视频。
该部分代码实现功能：循环将两幅图像拼接为全景图片，储存为有顺序的图像序列，方便后续拼成视频。
方法：以stit...
图像拼接(不投影到柱面)(渐入渐出融合) matlab程序
1，先拍摄一组图片，比如两幅图：A、B
我直接用网上的两幅图：
2，分别投影到柱面坐标系
就用自己写的柱面投影程序 matlab里 结果：
3，开始配准第一步：SIFT得到...
用matlab实现多张图片合并
I = imread('qiegray.jpg');
i=imrotate(I,45);
j=imrotate(I,315);
h1 = axes('position', [0.0 0.0 1.0 1...
基于matlab的图像拼接方法
图像拼接是一项应用广泛的图像处理技术。根据特征点的相互匹配，可以将多张小视角的图像拼接成为一张大视角的图像，在广角照片合成、卫星照片处理、医学图像处理等领域都有应用。早期的图像拼接主要...
matlab实现图像拼接
&& a=imread('512.bmp');
&& b=imread('731.bmp');
&& figure,imshow(a+b);
&& c=roipoly(a);//抠图，双击完成
...
matlab实现全景图像拼接技术
matlab实现全景图像拼接技术
Matlab将图像拼接成视频
%// Will open an avi file name test.avi in local folder
aviobj = avifile('test.avi');
%// the qual...
图像拼接算法总结（一）
图像的拼接技术包括三大部分：特征点提取与匹配、图像配准、图像融合。
1、基于SRUF 的特征点的提取与匹配
为了使拼接具有良好的精度和鲁棒性，同时又使其具有较好的实时性，本实验采用SURF 算法完成图...
没有更多推荐了，}