function imdst = process_corr(file, image1, image2)

M = load_corr(file);
% odriznou sloupce!!!
% M = M(:,8:end);

% evaluate transformation matrix
A = [];
for i = 1:size(M, 2),
    %B = [-M(1,i) -M(2,i) -1 0 0 0 M(1,i)*M(3,i) M(2,i)*M(3,i) M(3,i); ...
    %     0 0 0 -M(1,i) -M(2,i) -1 M(1,i)*M(4,i) M(2,i)*M(4,i) M(4,i)];
    B = [-M(3,i) -M(4,i) -1 0 0 0 M(3,i)*M(1,i) M(4,i)*M(1,i) M(1,i); ...
         0 0 0 -M(3,i) -M(4,i) -1 M(3,i)*M(2,i) M(4,i)*M(2,i) M(2,i)];
    A = [A; B];
end

[U,S,V] = svd(A);
h = V(:,end);

H = reshape(h, 3, 3)';

% process images
im1 = imread(image1);
im2 = imread(image2);
imdst = zeros(2*size(im1,1), 2*size(im1,2), 3);

xoffset = ceil(size(imdst,1)/4):ceil(3*size(imdst,1)/4) - 1;
yoffset = 1:ceil(size(im1,2));

imdst(xoffset, yoffset, :) = im1(:,:,:);

if 0
    xoff = ceil(size(imdst,1)/4);
	for i = 1:size(im2,1)
       	for j = 1:size(im2,2)
   	       G = H*[j;i;1];
           G = G/G(3,1);
           %ceil(size(imdst,1)/4) + round(G(2))
           %round(G(1))
           imdst(round(G(2)) + xoff, round(G(1)),:) = im2(i,j,:);
        end
    end
else
    xoff = ceil(size(imdst,1)/4);
    cr1 = H*[0;0;1];
    cr1 = cr1/cr1(3);
    cr2 = H*[size(im2,2);0;1];
    cr2 = cr2/cr2(3);
    cr3 = H*[0;size(im2,1);1];
    cr3 = cr3/cr3(3);
    cr4 = H*[size(im2,2);size(im2,1);1];
    cr4 = cr4/cr4(3);
    y1 = floor(min(min(cr1(1),cr2(1)), min(cr3(1),cr4(1))));
    y2 = ceil(max(max(cr1(1),cr2(1)), max(cr3(1),cr4(1))));
    x1 = floor(min(min(cr1(2),cr2(2)), min(cr3(2),cr4(2))));
    x2 = ceil(max(max(cr1(2),cr2(2)), max(cr3(2),cr4(2))));
	iH = inv(H);
    maxx = size(im2,1);
    maxy = size(im2,2);
    if 0
        for i = x1:x2
            for j = y1:y2
                G = iH*[j;i;1];
                G = G/G(3);
                x = round(G(2));
                y = round(G(1));
                if (x > 0 & x <= maxx & y > 20 & y <= maxy)
                    imdst(i + xoff, j, :) = im2(x, y, :);
          	    end
   	        end
        end
    else
        % bilinear interpolation
        for i = x1:x2
            for j = y1:y2
                G = iH*[j;i;1];
                G = G/G(3);
                x = floor(G(2));
                y = floor(G(1));
                if (x > 1 & x < maxx & y > 20 & y < maxy)
                    dx = G(2) - x;
                    dy = G(1) - y;
                    imdst(i + xoff, j, :) = (1 - dy)*((1 - dx)*double(im2(x, y, :)) + dx*double(im2(x + 1, y, :))) + dy*((1 - dx)*double(im2(x, y + 1, :)) + dx*double(im2(x + 1, y + 1, :)));
                    %imdst(i + xoff, j, 1) = (1 - dy)*((1 - dx)*im2(x, y, 1) + dx*im2(x + 1, y, 1)) + dy*((1 - dx)*im2(x, y + 1, 1) + dx*im2(x + 1, y + 1, 1));
                    %imdst(i + xoff, j, 2) = (1 - dy)*((1 - dx)*im2(x, y, 2) + dx*im2(x + 1, y, 2)) + dy*((1 - dx)*im2(x, y + 1, 2) + dx*im2(x + 1, y + 1, 2));
                    %imdst(i + xoff, j, 3) = (1 - dy)*((1 - dx)*im2(x, y, 3) + dx*im2(x + 1, y, 3)) + dy*((1 - dx)*im2(x, y + 1, 3) + dx*im2(x + 1, y + 1, 3));
          	    end
   	        end
        end
	end
end

imdst = imdst/255;
