パノラマの作成

ライブスクリプトを開く

この例では、複数のイメージを自動的に繋ぎ合わせてパノラマを作成する方法を示します。イメージを繋ぎ合わせる手続きは、特徴ベースのイメージレジストレーションを拡張したものです。1 組のイメージペアのレジストレーションを行う代わりに、複数のイメージペアを互いに比較しながらレジストレーションを連続的に行い、パノラマを形成します。

イメージの読み込み

この例で使用するイメージセットには建物の写真が含まれています。これらの写真は、キャリブレーションされていないスマートフォンカメラを左から右へ水平に移動させて、建物のすべての部分をカバーするように撮影したものです。イメージにはレンズ歪みの影響が比較的少ないため、カメラのキャリブレーションは必要ありません。しかし、レンズ歪みがある場合は、パノラマを作成する前にカメラのキャリブレーションを行い、イメージの歪みを補正しなければなりません。必要に応じてカメラキャリブレーターアプリを使用し、カメラのキャリブレーションを行います。

イメージを読み込んで表示します。

buildingDir = fullfile(toolboxdir("vision"),"visiondata","building");
buildingScene = imageDatastore(buildingDir);

montage(buildingScene.Files)

Figure contains an axes object. The hidden axes object contains an object of type image.

イメージペアのレジストレーション

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

$I (n)$ と $I (n - 1)$ の間で特徴を検出し、マッチングします。
$I (n)$ を $I (n - 1)$ にマッピングする幾何学的変換 $T (n)$ を推定します。
$I (n)$ をパノラマイメージに $T (1) * T (2) * . . . * T (n - 1) * T (n)$ としてマッピングする変換を計算します。

イメージセットから最初のイメージを読み取ります。

I = readimage(buildingScene,1);

最初のイメージで特徴を初期化します。

grayImage = im2gray(I);
points = detectSURFFeatures(grayImage);
[features,points] = extractFeatures(grayImage,points);

すべての変換を単位行列に初期化します。建物のイメージがカメラにかなり近いため、射影変換を使用していることに注意してください。より遠くから撮影されたシーンの場合は、アフィン変換の使用を検討してください。

numImages = numel(buildingScene.Files);
tforms(numImages) = projtform2d;

イメージサイズを保持する変数を初期化します。

imageSize = zeros(numImages,2);

残りのイメージペアについても同じ処理を繰り返します。各イメージについて、前のイメージのポイントと特徴を保存します。その後、それを読み取ってグレースケールに変換し、イメージサイズを保存します。イメージから SURF 特徴量を検出して抽出し、前のイメージとの対応関係を見つけます。次に、2 つの対応間の変換を推定します。

for n = 2:numImages
    pointsPrevious = points;
    featuresPrevious = features;
        
    I = readimage(buildingScene, n);
    grayImage = im2gray(I);    
    imageSize(n,:) = size(grayImage);
    
    points = detectSURFFeatures(grayImage);
    [features,points] = extractFeatures(grayImage,points);
  
    indexPairs = matchFeatures(features,featuresPrevious,Unique=true);
    matchedPoints = points(indexPairs(:,1), :);
    matchedPointsPrev = pointsPrevious(indexPairs(:,2), :);        
    
    tforms(n) = estgeotform2d(matchedPoints, matchedPointsPrev,...
        "projective",Confidence=99.9,MaxNumTrials=2000);
    
    tforms(n).A = tforms(n-1).A * tforms(n).A; 
end

この段階では、tforms 内のすべての変換は最初のイメージに対する相対的な変換になっています。これによりイメージの遂次処理が可能になり、イメージのレジストレーション手順のコーディングが簡単になります。ただし、最初のイメージからパノラマを開始すると、パノラマを構成するイメージのほとんどが歪む傾向があるため、見た目があまり美しくない結果になることがよくあります。シーンの中央で歪みが一番少なくなるように変換を調整することで、より見栄えの良いパノラマを作成することができます。この改善を行うには、中央のイメージの変換を反転させて、この逆変換をそれ以外のすべてのイメージに適用します。まずprojtform2d outputLimits メソッドを使用して各変換の出力範囲を求めます。その後、この出力範囲を使ってシーンのほぼ中央にあるイメージを自動的に見つけます。

各変換の出力範囲を計算します。

for idx = 1:numel(tforms)           
    [xlim(idx,:),ylim(idx,:)] = outputLimits(tforms(idx),[1 imageSize(idx,2)],[1 imageSize(idx,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 idx = 1:numel(tforms)    
    tforms(idx).A = Tinv.A * tforms(idx).A;
end

パノラマの初期化

すべてのイメージをマッピングする対象となる、初期の空白のパノラマを作成します。outputLimits メソッドを使用して、すべての変換にわたる出力範囲の最小値と最大値を計算します。これらの値は、パノラマのサイズを自動計算するために使用されます。

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

出力範囲の最小値と最大値を見つけます。

xMin = min([1; xlim(:)]);
xMax = max([maxImageSize(2); xlim(:)]);

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

パノラマの幅と高さを計算します。

width  = round(xMax - xMin);
height = round(yMax - yMin);

イメージ I のタイプと特性に合わせて、次元が [height width 3] であるゼロの配列を設定し、パノラマを空白のキャンバスで初期化します。

panorama = zeros([height width 3],"like",I);