このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

特徴に基づくパノラマ イメージの繋ぎ合わせ

この例では、特徴ベースのイメージ レジストレーションの手法を使用してパノラマを自動作成する方法を説明します。

概要

特徴の検出とマッチングは、イメージ レジストレーション、追跡、オブジェクト検出といった数多くのコンピューター ビジョン アプリケーションで使用される強力な手法です。この例では、特徴ベースの手法を使用して一連のイメージを自動的に繋ぎ合わせます。イメージを繋ぎ合わせる手続きは、特徴ベースのイメージ レジストレーションを拡張したものです。1 組のイメージ ペアのレジストレーションを行う代わりに、複数のイメージ ペアを互いに比較しながらレジストレーションを連続的に行い、パノラマを形成します。

手順 1 - イメージの読み込み

この例で使用するイメージ セットには建物の写真が含まれています。これらの写真は、キャリブレーションされていないスマート フォン カメラを左から右へ水平に移動させて、建物のすべての部分をカバーするように撮影したものです。

下に示すように、イメージにはレンズ歪みの影響が比較的少ないため、カメラのキャリブレーションは必要ありません。しかし、レンズ歪みがある場合は、パノラマを作成する前にカメラのキャリブレーションを行い、イメージの歪みを補正しなければなりません。必要に応じてカメラ キャリブレーターアプリを使用し、カメラのキャリブレーションを行うことができます。

% Load images.
buildingDir = fullfile(toolboxdir('vision'),'visiondata','building');
buildingScene = imageDatastore(buildingDir);

% Display images to be stitched.
montage(buildingScene.Files)

手順 2 - イメージ ペアのレジストレーション

パノラマを作成するには、まず次の手順によって、連続したイメージ ペアのレジストレーションを行います。

1. $$I(n)$$ と $$I(n-1)$$ の間で特徴を検出し、マッチングします。

2. $$I(n)$$ を $$I(n-1)$$ にマッピングする幾何学的変換 $$T(n)$$ を推定します。

3. $$I(n)$$ をパノラマ イメージに $$T(1)*T(2)*...*T(n-1)*T(n)$$ としてマッピングする変換を計算します。

% Read the first image from the image set.
I = readimage(buildingScene,1);

% Initialize features for I(1)
grayImage = im2gray(I);
points = detectSURFFeatures(grayImage);
[features, points] = extractFeatures(grayImage,points);

% Initialize all the transformations to the identity matrix. Note that the
% projective transformation is used here because the building images are fairly
% close to the camera. For scenes captured from a further distance, you can use
% affine transformations.
numImages = numel(buildingScene.Files);
tforms(numImages) = projtform2d;

% Initialize variable to hold image sizes.
imageSize = zeros(numImages,2);

% Iterate over remaining image pairs
for n = 2:numImages
    % Store points and features for I(n-1).
    pointsPrevious = points;
    featuresPrevious = features;
        
    % Read I(n).
    I = readimage(buildingScene, n);
    
    % Convert image to grayscale.
    grayImage = im2gray(I);    
    
    % Save image size.
    imageSize(n,:) = size(grayImage);
    
    % Detect and extract SURF features for I(n).
    points = detectSURFFeatures(grayImage);    
    [features, points] = extractFeatures(grayImage, points);
  
    % Find correspondences between I(n) and I(n-1).
    indexPairs = matchFeatures(features, featuresPrevious, 'Unique', true);
       
    matchedPoints = points(indexPairs(:,1), :);
    matchedPointsPrev = pointsPrevious(indexPairs(:,2), :);        
    
    % Estimate the transformation between I(n) and I(n-1).
    tforms(n) = estgeotform2d(matchedPoints, matchedPointsPrev,...
        'projective', 'Confidence', 99.9, 'MaxNumTrials', 2000);
    
    % Compute T(1) * T(2) * ... * T(n-1) * T(n).
    tforms(n).A = tforms(n-1).A * tforms(n).A; 
end

この時点では、tforms のすべての変換が最初のイメージを基準として行われます。これは、すべてのイメージを順次処理できるので、イメージ レジストレーションの手続きをコード化するには便利な方法です。ただし、最初のイメージをパノラマの始点として使用すると、パノラマを構成するイメージの大部分で歪みが生じる傾向があるため、最高に美しいパノラマは生成されません。シーンの中央で歪みが一番少なくなるよう変換を修正することで、見栄えの良いパノラマを作成することができます。これを行うには、中央のイメージの変換を反転させて、その変換をそれ以外のすべてのイメージに適用します。

まずprojtform2d outputLimits メソッドを使用して各変換の出力範囲を求めます。その後、この出力範囲を使ってシーンのほぼ中央にあるイメージを自動的に見つけます。

% Compute the output limits for each transformation.
for i = 1:numel(tforms)           
    [xlim(i,:), ylim(i,:)] = outputLimits(tforms(i), [1 imageSize(i,2)], [1 imageSize(i,1)]);    
end

次に、各変換について X の範囲の平均を求め、中央にあるイメージを見つけます。ここではシーンが横長であることがわかっているので、X の範囲のみを使用します。他のイメージ セットを使用する場合は、必要に応じて X と Y の両方の範囲を使用して、中央にあるイメージを見つけます。

avgXLim = mean(xlim, 2);
[~,idx] = sort(avgXLim);
centerIdx = floor((numel(tforms)+1)/2);
centerImageIdx = idx(centerIdx);

最後に、中央にあるイメージの逆変換を他のすべてのイメージに適用します。

Tinv = invert(tforms(centerImageIdx));
for i = 1:numel(tforms)    
    tforms(i).A = Tinv.A * tforms(i).A;
end

手順 3 - パノラマの初期化

次に、すべてのイメージをマッピングする対象となる、初期の空白のパノラマを作成します。

outputLimits メソッドを使用して、すべての変換にわたる出力範囲の最小値と最大値を計算します。これらの値は、パノラマのサイズを自動計算するために使用されます。

for i = 1:numel(tforms)           
    [xlim(i,:), ylim(i,:)] = outputLimits(tforms(i), [1 imageSize(i,2)], [1 imageSize(i,1)]);
end

maxImageSize = max(imageSize);

% Find the minimum and maximum output limits. 
xMin = min([1; xlim(:)]);
xMax = max([maxImageSize(2); xlim(:)]);

yMin = min([1; ylim(:)]);
yMax = max([maxImageSize(1); ylim(:)]);

% Width and height of panorama.
width  = round(xMax - xMin);
height = round(yMax - yMin);

% Initialize the "empty" panorama.
panorama = zeros([height width 3], 'like', I);

手順 4 - パノラマの作成

imwarpを使用してイメージをパノラマにマッピングし、vision.AlphaBlenderを使用してイメージを重ね合わせます。

blender = vision.AlphaBlender('Operation', 'Binary mask', ...
    'MaskSource', 'Input port');  

% Create a 2-D spatial reference object defining the size of the panorama.
xLimits = [xMin xMax];
yLimits = [yMin yMax];
panoramaView = imref2d([height width], xLimits, yLimits);

% Create the panorama.
for i = 1:numImages
    
    I = readimage(buildingScene, i);   
   
    % Transform I into the panorama.
    warpedImage = imwarp(I, tforms(i), 'OutputView', panoramaView);
                  
    % Generate a binary mask.    
    mask = imwarp(true(size(I,1),size(I,2)), tforms(i), 'OutputView', panoramaView);
    
    % Overlay the warpedImage onto the panorama.
    panorama = step(blender, panorama, warpedImage, mask);
end

figure
imshow(panorama)

まとめ

この例では、特徴ベースのイメージ レジストレーションの手法を使用してパノラマを自動作成する方法を説明しました。この例に追加の手法を取り入れて、パノラマ イメージのブレンドと配置を改善することもできます [1]。

参考文献

[1] Matthew Brown and David G. Lowe. 2007. Automatic Panoramic Image Stitching using Invariant Features.Int. J. Comput. Vision 74, 1 (August 2007), 59-73.