function RegImgSequence
% RegImgSequence register an image sequence by minimum spanning tree method

input_folder = 'Images\';
output_folder = 'Images_reg\';


files = dir([input_folder '*jpg']);
I = double(rgb2gray(imread([input_folder files(1).name])));
J = I;
N = length(files);

g_width = 5; % graph width

Wt = zeros(N,N); % edge weights of a graph
for i=1:N,
    I = double(rgb2gray(imread([input_folder files(i).name])));
    I(:) = medfilt2(I,[3 3]);
    for j=max(1,i-g_width):i-1,
        J(:) = double(rgb2gray(imread([input_folder files(j).name])));
        J(:) = medfilt2(J,[3 3]);
        Wt(i,j) = mean2(abs(imresize(I-J,0.25)));
        Wt(j,i) = Wt(i,j);
    end
    disp(['Building graph for image ' num2str(i) ' of ' num2str(N)]);
end
disp(' ');
% call minimum spanning tree
[~, pred]= graphminspantree(sparse(Wt));
% MST construction
Tr = zeros(N,N);
for i=1:N,
    if pred(i) ~= 0,
        Tr(i,pred(i)) = 1;
        Tr(pred(i),i) = 1;
    end
end

% decide on an anchor image
path = graphpred2path(pred, N);
anchor = path(floor(length(path)/2));
disp(['Anchor image: ' files(anchor).name]);
disp(' ');
% traverse the tree to register images; 
[disc, Tpred,~] = graphtraverse(sparse(Tr),anchor,'Method','BFS','Directed',false);
% always register a child to its parent
I = imread([input_folder files(anchor).name]);
J = I;
imwrite(J,[output_folder files(anchor).name],'jpg');
xform = eye(3);
trans = zeros(3,3,N);
trans(:,:,anchor) = xform;
for i=2:N,
    child = disc(i);
    I(:) = imread([input_folder files(child).name]);
    parent = Tpred(disc(i));
    J(:) = imread([input_folder files(parent).name]);
    xform = Reg2Im(I,J);
    xform(:) = xform*trans(:,:,parent);
    trans(:,:,child) = xform;
    tform = affine2d(xform);
    K = double(I);
    K(:,:,1) = imwarp(double(I(:,:,1)),tform,'OutputView',imref2d(size(I(:,:,1))));
    K(:,:,2) = imwarp(double(I(:,:,2)),tform,'OutputView',imref2d(size(I(:,:,2))));
    K(:,:,3) = imwarp(double(I(:,:,3)),tform,'OutputView',imref2d(size(I(:,:,3))));
    
    imwrite(K/max(K(:)),[output_folder files(child).name],'jpg');
    disp(['Registering image ' num2str(i) ' of ' num2str(N)]);
    
end

end

function xform = Reg2Im(J,I)
    % I: fixed image
    % J: moving image
    
    winsiz = 25;
    A = medfilt2(double(rgb2gray(J)),[5 5]);
    A = A(winsiz+1:end-winsiz,winsiz+1:end-winsiz);
    B = medfilt2(double(rgb2gray(I)),[5 5]);
    c = normxcorr2(A,B);
    [ypeak, xpeak] = find(c(2*winsiz:end-2*winsiz,2*winsiz:end-2*winsiz)==max(max(c(2*winsiz:end-2*winsiz,2*winsiz:end-2*winsiz))));
    yoffSet = ypeak(1) - size(A,1) + winsiz - 1;
    xoffSet = xpeak(1) - size(A,2) + winsiz - 1;
    
    xform=eye(3);
    xform(3,1) = xoffSet;
    xform(3,2) = yoffSet;    
end
