% ortofoto demo - data upravena z klad55.jpg 
% http://geoportal.cuzk.cz/UKAZKOVA_DATA/ORTOFOTO_JPG.zip
% ukazka primeho vypoctu projektivni trasformace

url = 'http://geo.fsv.cvut.cz/user/gin/git1/matlab/img/';
A = imread(strcat(url, 'klad55-ortofoto.jpg')); % ortofoto snimek (vyrez)      
B = imread(strcat(url, 'klad55-perspektiva.jpg'));

A_ib = [115 123; 445 177; 467 302; 332 362; 
        279 446;  95 234; 389 465; 170 533];
B_ib = [ 42  93; 409  70; 467 174; 345 260; 
        312 353;  48 200; 436 337; 214 467];
%{    
B( 1 , : ,:) = 0;    B( 1 , : ,1) = 255;    
B(end, : ,:) = 0;    B(end, : ,1) = 255; 
B( : , 1 ,:) = 0;    B( : , 1 ,1) = 255;    
B( : ,end,:) = 0;    B( : ,end,1) = 255;
%}
     
M = zeros(size(A_ib, 1)*2, 8);   % rovnice oprav
r = zeros(size(M, 1),1);
j = 1;
for i = 1:size(A_ib,1)
    X = A_ib(i,1);   Y = A_ib(i,2);
    x = B_ib(i,1);   y = B_ib(i,2);
    M(j,1) = x;  M(j,2) = y;  M(j,3) = 1;  M(j,7) = -x*X;  M(j,8) = -y*X;     
    r(j)   = X;  j = j + 1;
    M(j,4) = x;  M(j,5) = y;  M(j,6) = 1;  M(j,7) = -x*Y;  M(j,8) = -y*Y;     
    r(j)   = Y;  j = j + 1;
end
scale = max(abs(M));   scale(eq(scale,0)) = 1;
%scale(:) = 1.0
for j = 1:size(M,2)
    M(:,j) = M(:,j)./scale(j); 
end;
%cond(M)                         % cislo podminenosti

[U, W, V] = svd(M, 'econ');      % reseni MNC
%q = V*inv(W)*U'*r;              % inverzi nahradime prevracenou hodnotou
w = diag(W);
w(eq(w,0)) = 1;                  % linearne zavisle sloupce
w = 1./w;
q = V*diag(w)*U'*r;              % mat nasobeni pouze pro strucnejsi zapis
a = q(1)/scale(1); b = q(2)/scale(2); c = q(3)/scale(3); 
d = q(4)/scale(4); e = q(5)/scale(5); f = q(6)/scale(6); 
g = q(7)/scale(7); h = q(8)/scale(8);
%[a b c d e f g h]               % parametry projektivni transformace

R = size(A,1);
S = size(A,2);
for r = 1:size(B,1)
    for s = 1:size(B,2)
       q =  g*s + h*r + 1;
       x = round( (a*s + b*r + c)/q );
       y = round( (d*s + e*r + f)/q );
       if x < 1 || x > S || y < 1 || y > R
          continue;
       end
       A(y,x,1) = B(r,s,1);
       A(y,x,2) = B(r,s,2);
       A(y,x,3) = B(r,s,3)*0;    % #### zvyrazneni vysledku ####
    end
end

figure(1), imshow(A);
imwrite(A, 'klad55-vystup-svd-m.png', 'png');