Zošívanie obrazu

1. Zvoľ sieľovú projektciu: rovná plocha
   
2. Označ obrázky ako \(I(1), ...I(5)\)
   
3. Pre každú dvojicu obrazu \(I(n), ...I(n-1)\) vykonaj registráciu, nájdi transformáciu \(I(n)\) do súradnicového systému \(I(n-1)\)
1. Pomocou metódy SURF identifikuj význačné body na  \(I(n)\) aj \(I(n-1)\) . MATLAB : detectSURFFeatures() + extractFeatures()
2. Hľadaj medzi nimi jedinečné páry s najlepšou zhodou (MATLAB: matchFeatures(), s voľbou "Unique") na základe metriky SAD (Sum of Absolute Differences)
3. Na základe párov bodov odhadni transformáciu  \(I(n)\) do súradnicového systému \(I(n-1)\) pomocou metódy  MSAC (M-estimator sample consensus) ktorá je variantom známej metódy RANSAC (random sample consensus), pričom sa hľadá projektívna transformácia (MATLAB: estgeotform2d() s parametrom "projective") 
4. vypočítaj výslednú transformáciu \(I(n)\) voči počiatočnému \(I(1)\) ako  \(\prod_{n=1}^n T(n)\) 
4. Ak je centrálny obraz iný ako \(I(1)\) prepočítaj\(T(n)\) voči tomuto obrazu
5. K \(I(1)\) postupne pre prišívaj ďalšie obrázky, pričom ako šev použi transformovanú hranu prišívaného obrazu a použi najjednoduchšie zmiešavanie t.j. priši rovno plné hodnoty pixelov v oblasti nového obrazu (MATLAB: vision.AlphaBlender, pričom operácia je "Binary mask")

Na rozdiel od situácie na Obr. 2, kde zošívanie prebiehalo horizontálne, v reálnych panorámach sa často zošíva aj v druhom smere, teda výsledný obraz je tvorený napr. z 20x10 vstupných prekrývajúcich sa obrazov. Na zošívanie obrazu existuje množstvo špecializovaných programov (výborný prehľad je uvedený napr. v [2])

Referencie

[1] Mathworks, Feature Based Panoramic Image Stitching, http://www.mathworks.com/help/vision/ug/feature-based-panoramic-image-stitching.html

[2] Comparison of photo stitching software, wikipedia, https://en.wikipedia.org/wiki/Comparison_of_photo_stitching_software#:~:text=Photo%20stitching%20software%20produce%20panoramic,more%20or%20less%20smooth%20image.

Prílohy

Zdrojový kód programu, ktorý vytvoril Obr. 2 aj Obr. 3

  clear all;close all;clc;
%https://www.mathworks.com/help/vision/ug/feature-based-panoramic-image-stitching.html
% Load images.
buildingDir = fullfile(toolboxdir('vision'),'visiondata','building');
buildingScene = imageDatastore(buildingDir);

[tforms,imageSize]=getTformsSizes(buildingScene);
panorama0=getPanorama(buildingScene,tforms,imageSize,0);
panorama1=getPanorama(buildingScene,tforms,imageSize,1);

fig=figure;
subplot(5,1,[1 1]);
montage(buildingScene.Files,"Size",[1 5], "BorderSize", [2 2],"BackgroundColor","w");
subplot(5,1,[2 3]);
imshow(panorama0)
subplot(5,1,[4 5]);
imshow(panorama1)
fig.PaperUnits = 'inches';
fig.PaperPosition = [0 0 10 12];
print(fig,"stitching.png",'-dpng');


%=========================================================
function [tforms,imageSize]=getTformsSizes(buildingScene)

    % 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
end

%=========================================================
function panorama = getPanorama(buildingScene,tforms,imageSize,doCentering)

    if(doCentering)
        % 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
        
        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
    end
    
    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.
    I = readimage(buildingScene,1);
    panorama = zeros([height width 3], 'like', I);
    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);
    numImages = numel(buildingScene.Files);
    % Create the panorama.
    for i = 1:numImages;
        
        I = readimage(buildingScene, i);   
        imwrite(I,sprintf("img%02i.png",0*5+i));    
        % Transform I into the panorama.
        warpedImage = imwarp(I, tforms(i), 'OutputView', panoramaView);
        imwrite(warpedImage,sprintf("img%02i.png",1*5+i));    
        % Generate a binary mask.    
        mask = imwarp(true(size(I,1),size(I,2)), tforms(i), 'OutputView', panoramaView);   
        imwrite(mask,sprintf("img%02i.png",2*5+i));
        % Overlay the warpedImage onto the panorama.
        panorama = step(blender, panorama, warpedImage, mask);
        imwrite(panorama,sprintf("img%02i.png",3*5+i));
    end
    fig=figure
    imds = imageDatastore(fullfile(".","img*"));
    montage(imds,"Size", [4 5],"BorderSize", [1 1],"BackgroundColor","w")
    print(fig,sprintf("stitchingDetailsDC%u.png",doCentering),'-dpng');
end